This notebook contains all multivariate analyses of zoobenthic community structure using the new, nearly unheard-of modeling methods: packages mvabund, boral.
Again, to make it self-contained, there will be the same repetitive setup/data import/preparation part.
Setup!
Define the working subdirectories.
## print the working directory, just to be on the safe side
paste("You are here: ", getwd())
data.dir <- "data" ## input data files
functions.dir <- "R" ## functions & scripts
save.dir <- "output" ## clean data, output from models & more complex calculations
figures.dir <- "figs" ## plots & figures
Import libraries.
library(here) ## painless relative paths to subdurectories, etc.
library(tidyverse) ## data manipulation, cleaning, aggregation
library(viridis) ## smart & pretty colour schemes
library(mvabund) ## multivariate modeling analyses in ecology
library(boral) ## more multivariate modeling analyses in ecology
Organize some commonly-used ggplot2 modifications into a more convenient (and less repetitive) format. One day, I MUST figure out the proper way to set the theme..
## ggplot settings & things that I keep reusing
# ggplot_theme <- list(
# theme_bw(),
# theme(element_text(family = "Times"))
# )
## always use black-and-white theme
theme_set(theme_bw())
## helper to adjust ggplot text size & avoid repetitions
text_size <- function(text.x = NULL,
text.y = NULL,
title.x = NULL,
title.y = NULL,
legend.text = NULL,
legend.title = NULL,
strip.x = NULL,
strip.y = NULL) {
theme(axis.text.x = element_text(size = text.x),
axis.text.y = element_text(size = text.y),
axis.title.x = element_text(size = title.x),
axis.title.y = element_text(size = title.y),
legend.text = element_text(size = legend.text),
legend.title = element_text(size = legend.title),
strip.text.x = element_text(size = strip.x),
strip.text.y = element_text(size = strip.y)
)
}
## log y/min + 1 transform - useful for species counts/biomass data visualization
log_y_min <- function(y) {
log(y / min(y[y > 0]) + 1)
}
Import zoobenthic abundance data (cleaned and prepared).
zoo.abnd.sand <- read_csv(here(save.dir, "abnd_sand_orig_clean.csv"))
## convert station to factor (better safe than sorry later, when the stations are not plotted in the order I want them)
(zoo.abnd.sand <- zoo.abnd.sand %>%
mutate(station = factor(station, levels = c("Kraimorie", "Chukalya", "Akin", "Sozopol", "Agalina", "Paraskeva")))
)
Remove the all-0 species (= not present in the current dataset).
Maybe also remove the singletons (species appearing only once in the whole dataset and represented by a single individual = so rare that it’s unlikely they carry important information, but it would probably improve the run times).
(zoo.abnd.flt.sand <- zoo.abnd.sand %>%
select(-c(station:replicate)) %>%
select(which(colSums(.) > 0))
)
Perform a model-based unconstrained ordination by fiting a pure latent variable model (package boral - Hui et al., 2014). This will allow to visualize the multivariate stations x species data - similar to nMDS, can be interpreted in the same way.
I’m including a (fixed) row effect to account for differences in site total abundance - this way, the ordination is in terms of species composition.
NB this takes about a million years to run!
lvm.sand <- boral(y = zoo.abnd.flt.sand,
family = "negative.binomial",
## we want to control for site effects - there are 6 sites with 9 replicates each
row.eff = "fixed", row.ids = matrix(rep(1:6, each = 9), ncol = 1),
## 2 latent variables = 2 axes on which to represent the zoobenthic data
lv.control = list(num.lv = 2)
# ## example control structure, to check if function does what I want, because otherwise it takes an intolerably long time, and I'll shoot myself if I have to wait for it again
# mcmc.control = list(n.burnin = 10, n.iteration = 100,
# n.thin = 1)
#
)
Check the summary and diagnostic plots for the LVM.
plot(lvm.sand)
NULL
The residuals plots look fine (no patterns in the residuals vs fitted, so variance is homogeneous, the quantile plot shows a normal distribution of the residuals) - the model fits the data pretty well.
Save the sand LVM.
write_rds(lvm.sand,
here(save.dir, "lvm_sand.RDS"))
Examine the biplot obtained by fitting the LVM, as well as the 20 most “important” species.
lvsplot(lvm.sand, jitter = T, biplot = TRUE, ind.spp = 20)
Only the first 20 ``most important'' latent variable coefficients included in biplot
All in all, the final result resembles the nMDS ordination very much - same 4 clusters (Kraimorie + Chukalya, AKin, Agalina, Sozopol + Paraskeva). Kraimorie and Chukalya are better distinguished on the LVM plot than on the MDS, but still.
The run time is extremely, extremely long (~1h), but the data don’t need to be transformed, and the model fit can be examined and adjusted if necessary.
The species singled out as significant are probably somewhat different - have to check!
Redo the biplot, because this one is not very pretty. I’m not adding the species on top, first because I’m too lazy to figure out the procedure for ordering them, and second because the plot gets too busy.
Make the plot and save it.
## save the LVM plot for the sand stations
ggsave(file = here(figures.dir, "lvm_sand.png"),
plot.lvm.sand,
width = 15, units = "cm", dpi = 300)
Let’s fit GLMs to the sites x species matrix to try and explain the observed differences in community structure by the variation of the environmental parameters.
These functions all come from package mvabund.
Import the environmental data - the one cleaned, prepared and saved in the previous notebook (classical multivariate methods). It contains long-term averages for the water column data (2009-2011 + 2013-2014) at each station, repeated for each replicate, and the sediment data (2013-2014), again repeated to the same number of replicates. Only the variables determined to be significant by PCA are kept.
env.sand <- read_csv(here(save.dir, "env_data_ordinations_sand.csv"))
## convert station to factor
(env.sand <- env.sand %>%
mutate(station = factor(station,
levels = c("Kraimorie", "Chukalya", "Akin", "Sozopol", "Agalina", "Paraskeva")))
)
Station is a factor, the rest of the variables are numeric.
Turn the zoobenthic data (minus the all-0 taxa) into a matrix - easier for the mvabund package and methods to deal with.
## there is already one subset of filtered count data (54 x 147) - use it
zoo.mvabnd.sand <- mvabund(zoo.abnd.flt.sand)
First, let’s see if the groups from the latent variable model (more or less equal to the clusters from the classical ordination) are valid, and which species exhibit a response.
## construct the vector of the clusters by hand, it's easier that way..
lvm.clusters.sand <- c(rep(1, times = 18), rep(2:4, each = 9), rep(3, times = 9))
## convert to factor
(lvm.clusters.sand <- factor(lvm.clusters.sand))
Check the model assumptions. 1. Mean-variance assumption => determines the choice of family parameter. Can be checked by plotting residuals vs fits: if little pattern - the chosen mean-variance assumption is plausible.
Another way: direct plotting (variance ~ mean), for each species within each factor level.
plot(manyglm(zoo.mvabnd.sand ~ lvm.clusters.sand, family = "negative.binomial"))
meanvar.plot(zoo.mvabnd.sand ~ lvm.clusters.sand, table = TRUE)
START SECTION 2
Plotting if overlay is TRUE
using grouping variable lvm.clusters.sand 290 mean values were 0 and could
not be included in the log-plot
using grouping variable lvm.clusters.sand 290 variance values were 0 and could not
be included in the log-plot
FINISHED SECTION 2
$mean
Chamelea.gallina Protodorvillea.kefersteini Oligochaeta Microdeutopus.versiculatus
1 5.0555556 22.7222222 29.2222222 0.1666667
2 10.1111111 190.5555556 187.0000000 0.0000000
3 140.3888889 0.8333333 0.2222222 0.0000000
4 0.3333333 60.2222222 20.2222222 136.5555556
Heteromastus.filiformis Melinna.palmata Prionospio.cirrifera Lentidium.mediterraneum
1 62.72222222 51.55555556 26.0000000 0.00000
2 0.00000000 0.00000000 3.0000000 0.00000
3 0.05555556 0.05555556 0.3333333 21.77778
4 0.22222222 0.00000000 0.0000000 0.00000
Eurydice.dollfusi Polygordius.neapolitanus Branchiostoma.lanceolatum Ampelisca.diadema
1 0.0000000 1.333333 0.0000000 9.944444
2 6.5555556 8.444444 23.1111111 6.666667
3 0.7777778 0.000000 0.5555556 2.111111
4 34.7777778 22.888889 9.5555556 1.333333
Melita.palmata Bittium.reticulatum Capitella.minima Nemertea Micronephthys.stammeri
1 0.00000 9.7222222 6.4444444 2.222222 3.944444
2 0.00000 3.8888889 0.1111111 3.555556 0.000000
3 0.00000 1.0555556 2.6666667 2.055556 3.277778
4 31.77778 0.8888889 2.1111111 6.666667 3.222222
Pseudocuma.longicorne Polycirrus.caliendrum Sphaerosyllis.hystrix
1 0.0000000 0.00000 1.0000000
2 3.3333333 0.00000 2.7777778
3 7.0000000 0.00000 0.2222222
4 0.1111111 17.22222 11.1111111
Monocorophium.acherusicum Polycirrus.jubatus Ophelia.limacina Spio.filicornis
1 4.9444444 0.00000 0.0000000 0.05555556
2 0.7777778 0.00000 12.3333333 0.00000000
3 0.1111111 0.00000 0.1111111 7.11111111
4 5.3333333 14.88889 2.2222222 0.00000000
Hirudinea Lindrilus.flavocapitatus Polydora.ciliata Streptosyllis.bidentata
1 0.05555556 0.2222222 6.7222222 0.00000
2 4.11111111 13.6666667 0.0000000 0.00000
3 0.00000000 0.0000000 0.1111111 0.00000
4 9.88888889 0.0000000 0.0000000 13.66667
Acari Anadara.kagoshimensis Abra.alba Parvicardium.exiguum Exogone.naidina
1 0.05555556 5.666667 3.0000000 4.0000000 1.2777778
2 0.33333333 0.000000 0.4444444 0.3333333 0.6666667
3 0.00000000 0.000000 2.3888889 0.9444444 1.2222222
4 11.88888889 0.000000 0.0000000 0.5555556 4.5555556
Nototropis.guttatus Microphthalmus.fragilis Lagis.koreni Nephtys.cirrosa
1 0.4444444 0.000000 4.2222222 1.0000000
2 4.7777778 6.333333 0.1111111 0.0000000
3 0.6111111 0.000000 0.0000000 3.1111111
4 2.3333333 2.777778 0.2222222 0.2222222
Lucinella.divaricata Perioculodes.longimanus Bathyporeia.guilliamsoniana
1 0.000000 2.1666667 0.0000000
2 0.000000 1.3333333 4.7777778
3 3.666667 0.5555556 1.1111111
4 0.000000 0.5555556 0.1111111
Glycera.tridactyla Tellina.tenuis Aonides.paucibranchiata Caecum.armoricum
1 2.7222222 0.5000000 2.8333333 0.0000000
2 0.0000000 0.5555556 0.0000000 0.2222222
3 0.7777778 2.3333333 0.0000000 0.0000000
4 0.1111111 0.0000000 0.4444444 5.7777778
Bodotria.arenosa Spisula.subtruncata Harmothoe.reticulata Eunice.vittata Turbellaria
1 0.3333333 0.05555556 1.888889 0.1666667 1.1111111
2 1.4444444 0.00000000 0.000000 0.0000000 0.4444444
3 0.0000000 2.72222222 0.000000 0.0000000 0.3888889
4 3.5555556 0.00000000 1.555556 4.4444444 1.0000000
Diogenes.pugilator Loripes.orbiculatus Mytilaster.lineatus Genetyllis.tuberculata
1 0.9444444 0.05555556 1.0555556 0.1111111
2 0.7777778 0.00000000 1.0000000 0.0000000
3 0.5555556 2.00000000 0.3333333 0.5000000
4 0.4444444 0.00000000 0.3333333 2.7777778
Iphinoe.tenella Dinophilus.gyrociliatus Perinereis.cultrifera Pitar.rudis
1 1.888889 0.27777778 0.05555556 1.2222222
2 0.000000 0.00000000 0.33333333 0.0000000
3 0.000000 0.05555556 0.27777778 0.0000000
4 0.000000 3.00000000 2.33333333 0.7777778
Syllis.hyalina Leiochone.leiopygos Salvatoria.clavata Amphibalanus.improvisus
1 0.000000 1.4444444 0.0000000 0.9444444
2 0.000000 0.0000000 0.5555556 0.4444444
3 0.000000 0.0000000 0.0000000 0.2222222
4 3.222222 0.2222222 2.5555556 0.2222222
Capitella.capitata Syllis.gracilis Decapoda.larvae Magelona.papillicornis
1 0.7777778 0.3888889 1.2777778 0.0000000
2 0.2222222 0.0000000 0.1111111 0.1111111
3 0.3888889 0.0000000 0.0000000 1.2222222
4 0.4444444 2.1111111 0.1111111 0.0000000
Tritia.neritea Alitta.succinea Eteone.sp. Schistomeringos.rudolphi Microphthalmus.sp.
1 0.000000 0.50000000 0.2222222 0.00000000 0.0000000
2 0.000000 0.11111111 0.0000000 0.00000000 0.0000000
3 1.222222 0.05555556 0.2222222 0.05555556 0.3888889
4 0.000000 1.11111111 1.4444444 1.88888889 1.0000000
Caecum.trachea Donax.venustus Mysta.picta Glycera.sp. Odostomia.plicata
1 0.0000000 0.0000000 0.1111111 0.6111111 0.1111111
2 0.5555556 0.0000000 0.2222222 0.0000000 0.0000000
3 0.0000000 0.8333333 0.4444444 0.0000000 0.6111111
4 1.1111111 0.0000000 0.3333333 0.3333333 0.0000000
Apseudopsis.ostroumovi Cytharella.costulata Kellia.suborbicularis
1 0.38888889 0.61111111 0.6666667
2 0.00000000 0.00000000 0.0000000
3 0.05555556 0.05555556 0.0000000
4 0.44444444 0.00000000 0.0000000
Leptosynapta.inhaerens Spionidae Tritia.reticulata Gastrosaccus.sanctus
1 0.000000 0.3888889 0.3888889 0.05555556
2 0.000000 0.1111111 0.5555556 0.66666667
3 0.000000 0.0000000 0.0000000 0.22222222
4 1.333333 0.4444444 0.0000000 0.00000000
Micromaldane.ornithochaeta Phoronida Holothuroidea Nereis.perivisceralis
1 0.0000000 0.55555556 0.000000 0
2 0.0000000 0.00000000 0.000000 0
3 0.6111111 0.05555556 0.000000 0
4 0.0000000 0.00000000 1.111111 1
Phyllodoce.sp. Amphitritides.gracilis Gammaridae Amphipoda Harmothoe.imbricata
1 0.11111111 0.0000000 0.00000000 0.0000000 0.3888889
2 0.22222222 0.0000000 0.77777778 0.7777778 0.0000000
3 0.05555556 0.0000000 0.05555556 0.0000000 0.0000000
4 0.44444444 0.8888889 0.00000000 0.0000000 0.0000000
Mytilus.galloprovincialis Pholoe.inornata Pisione.remota Ampithoe.sp.
1 0.05555556 0.0000000 0.05555556 0.0000000
2 0.00000000 0.0000000 0.44444444 0.3333333
3 0.00000000 0.0000000 0.00000000 0.1666667
4 0.66666667 0.7777778 0.22222222 0.0000000
Paradoneis.harpagonea Platyhelminthes Platynereis.dumerilii Rissoa.membranacea
1 0.0000000 0.0000000 0.0000000 0.0000000
2 0.0000000 0.2222222 0.0000000 0.0000000
3 0.3333333 0.0000000 0.0000000 0.3333333
4 0.0000000 0.4444444 0.6666667 0.0000000
Brachynotus.sexdentatus Eulalia.viridis Microdeutopus.gryllotalpa Nereis.pelagica
1 0.2777778 0.05555556 0.1111111 0.0000000
2 0.0000000 0.00000000 0.1111111 0.0000000
3 0.0000000 0.00000000 0.1111111 0.0000000
4 0.0000000 0.44444444 0.0000000 0.5555556
Rapana.venosa Sabellaria.taurica Steromphala.divaricata Aricidea.claudiae
1 0.2777778 0.2777778 0.0000000 0.2222222
2 0.0000000 0.0000000 0.0000000 0.0000000
3 0.0000000 0.0000000 0.0000000 0.0000000
4 0.0000000 0.0000000 0.5555556 0.0000000
Monocorophium.insidiosum Nereis.pulsatoria Parthenina.interstincta Protodrilus.sp.
1 0.0000000 0.16666667 0.05555556 0.0000000
2 0.0000000 0.00000000 0.00000000 0.4444444
3 0.2222222 0.05555556 0.16666667 0.0000000
4 0.0000000 0.00000000 0.00000000 0.0000000
Sternaspis.scutata Tellina.fabula Cerastoderma.edule Chondrochelia.savignyi
1 0.0000000 0.0000000 0.0000000 0.0000000
2 0.4444444 0.0000000 0.0000000 0.0000000
3 0.0000000 0.2222222 0.1666667 0.0000000
4 0.0000000 0.0000000 0.0000000 0.3333333
Crassicorophium.crassicorne Magelona.rosea Polititapes.aureus Polychaeta.larvae
1 0.0000000 0.0000000 0.1666667 0.0000000
2 0.3333333 0.0000000 0.0000000 0.3333333
3 0.0000000 0.1666667 0.0000000 0.0000000
4 0.0000000 0.0000000 0.0000000 0.0000000
Retusa.variabilis Brachystomia.scalaris Eteone.flava Hydrobia.sp. Mactra.stultorum
1 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000
2 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000
3 0.1666667 0.1111111 0.1111111 0.05555556 0.1111111
4 0.0000000 0.0000000 0.0000000 0.11111111 0.0000000
Nephtyidae Papillicardium.papillosum Abra.sp. Acanthocardia.tuberculata Actiniaria
1 0.0000000 0.0000000 0.00000000 0.00000000 0.00000000
2 0.0000000 0.0000000 0.00000000 0.00000000 0.00000000
3 0.1111111 0.1111111 0.05555556 0.05555556 0.05555556
4 0.0000000 0.0000000 0.00000000 0.00000000 0.00000000
Cardiidae Cerastoderma.glaucum Cumopsis.goodsir Elaphognathia.bacescoi
1 0.00000000 0.05555556 0.00000000 0.0000000
2 0.00000000 0.00000000 0.00000000 0.0000000
3 0.05555556 0.00000000 0.05555556 0.0000000
4 0.00000000 0.00000000 0.00000000 0.1111111
Eurydice.pontica Hydrobia.acuta Iphinoe.sp. Lysidice.ninetta Maldanidae
1 0.0000000 0.05555556 0.0000000 0.0000000 0.00000000
2 0.1111111 0.00000000 0.1111111 0.1111111 0.00000000
3 0.0000000 0.00000000 0.0000000 0.0000000 0.05555556
4 0.0000000 0.00000000 0.0000000 0.0000000 0.00000000
Microdeutopus.sp. Neanthes.sp. Nereididae Pestarella.candida Pisces.larvae
1 0.00000000 0.05555556 0.0000000 0.0000000 0.05555556
2 0.00000000 0.00000000 0.1111111 0.0000000 0.00000000
3 0.05555556 0.00000000 0.0000000 0.0000000 0.00000000
4 0.00000000 0.00000000 0.0000000 0.1111111 0.00000000
Retusa.truncatula Rhithropanopeus.harrisii Stenothoe.monoculoides Upogebia.pusilla
1 0.00000000 0.0000000 0.0000000 0.0000000
2 0.00000000 0.0000000 0.0000000 0.0000000
3 0.05555556 0.0000000 0.0000000 0.0000000
4 0.00000000 0.1111111 0.1111111 0.1111111
$var
Chamelea.gallina Protodorvillea.kefersteini Oligochaeta Microdeutopus.versiculatus
1 32.05556 707.0359477 3.339477e+02 0.500
2 38.86111 7209.2777778 3.244850e+04 0.000
3 7120.36928 0.9705882 3.006536e-01 0.000
4 0.25000 2202.9444444 1.089444e+02 8235.778
Heteromastus.filiformis Melinna.palmata Prionospio.cirrifera Lentidium.mediterraneum
1 1.564683e+03 2.225908e+03 339.2941176 0.000
2 0.000000e+00 0.000000e+00 17.7500000 0.000
3 5.555556e-02 5.555556e-02 0.4705882 2097.242
4 4.444444e-01 0.000000e+00 0.0000000 0.000
Eurydice.dollfusi Polygordius.neapolitanus Branchiostoma.lanceolatum Ampelisca.diadema
1 0.000000 10.23529 0.0000000 107.93791
2 46.277778 50.02778 52.8611111 59.50000
3 3.712418 0.00000 0.7320261 7.51634
4 638.944444 68.11111 12.0277778 1.50000
Melita.palmata Bittium.reticulatum Capitella.minima Nemertea Micronephthys.stammeri
1 0.0000 142.447712 48.2614379 4.418301 12.173203
2 0.0000 9.861111 0.1111111 45.777778 0.000000
3 0.0000 1.349673 38.0000000 9.702614 5.388889
4 281.4444 1.611111 5.3611111 62.750000 3.944444
Pseudocuma.longicorne Polycirrus.caliendrum Sphaerosyllis.hystrix
1 0.0000000 0.0000 2.7058824
2 16.0000000 0.0000 8.9444444
3 38.0000000 0.0000 0.5359477
4 0.1111111 323.6944 72.8611111
Monocorophium.acherusicum Polycirrus.jubatus Ophelia.limacina Spio.filicornis
1 17.7026144 0.00000 0.0000000 0.05555556
2 1.4444444 0.00000 42.5000000 0.00000000
3 0.1045752 0.00000 0.1045752 136.81045752
4 26.0000000 64.11111 6.1944444 0.00000000
Hirudinea Lindrilus.flavocapitatus Polydora.ciliata Streptosyllis.bidentata
1 0.05555556 0.8888889 64.5653595 0.00
2 8.11111111 338.2500000 0.0000000 0.00
3 0.00000000 0.0000000 0.1045752 0.00
4 36.61111111 0.0000000 0.0000000 18.75
Acari Anadara.kagoshimensis Abra.alba Parvicardium.exiguum Exogone.naidina
1 0.05555556 10.23529 10.470588 7.0588235 1.506536
2 0.50000000 0.00000 1.777778 0.2500000 0.750000
3 0.00000000 0.00000 6.604575 1.2320261 1.947712
4 69.61111111 0.00000 0.000000 0.5277778 59.277778
Nototropis.guttatus Microphthalmus.fragilis Lagis.koreni Nephtys.cirrosa
1 0.3790850 0.00000 7.3594771 3.2941176
2 8.1944444 33.50000 0.1111111 0.0000000
3 0.8398693 0.00000 0.0000000 6.4575163
4 5.2500000 13.94444 0.1944444 0.4444444
Lucinella.divaricata Perioculodes.longimanus Bathyporeia.guilliamsoniana
1 0.00000 2.8529412 0.0000000
2 0.00000 3.2500000 6.9444444
3 56.23529 0.7320261 1.7516340
4 0.00000 1.7777778 0.1111111
Glycera.tridactyla Tellina.tenuis Aonides.paucibranchiata Caecum.armoricum
1 14.0947712 0.5000000 18.852941 0.0000000
2 0.0000000 0.5277778 0.000000 0.4444444
3 1.8300654 16.5882353 0.000000 0.0000000
4 0.1111111 0.0000000 1.027778 55.9444444
Bodotria.arenosa Spisula.subtruncata Harmothoe.reticulata Eunice.vittata Turbellaria
1 0.7058824 0.05555556 2.928105 0.1470588 4.3398693
2 5.2777778 0.00000000 0.000000 0.0000000 0.5277778
3 0.0000000 6.44771242 0.000000 0.0000000 0.7222222
4 4.5277778 0.00000000 1.277778 13.5277778 1.2500000
Diogenes.pugilator Loripes.orbiculatus Mytilaster.lineatus Genetyllis.tuberculata
1 0.7614379 0.05555556 0.7614379 0.1045752
2 1.1944444 0.00000000 2.0000000 0.0000000
3 0.3790850 9.05882353 0.2352941 0.8529412
4 0.2777778 0.00000000 0.2500000 15.1944444
Iphinoe.tenella Dinophilus.gyrociliatus Perinereis.cultrifera Pitar.rudis
1 3.51634 0.56535948 0.05555556 8.771242
2 0.00000 0.00000000 1.00000000 0.000000
3 0.00000 0.05555556 0.21241830 0.000000
4 0.00000 5.75000000 7.50000000 1.194444
Syllis.hyalina Leiochone.leiopygos Salvatoria.clavata Amphibalanus.improvisus
1 0.000000 2.0261438 0.0000000 0.8790850
2 0.000000 0.0000000 0.7777778 0.5277778
3 0.000000 0.0000000 0.0000000 0.1830065
4 8.444444 0.4444444 11.5277778 0.1944444
Capitella.capitata Syllis.gracilis Decapoda.larvae Magelona.papillicornis
1 1.5947712 0.4869281 14.8006536 0.0000000
2 0.4444444 0.0000000 0.1111111 0.1111111
3 0.3692810 0.0000000 0.0000000 0.8888889
4 0.7777778 2.8611111 0.1111111 0.0000000
Tritia.neritea Alitta.succinea Eteone.sp. Schistomeringos.rudolphi Microphthalmus.sp.
1 0.000000 0.85294118 0.3006536 0.00000000 0.0000000
2 0.000000 0.11111111 0.0000000 0.00000000 0.0000000
3 3.477124 0.05555556 0.3006536 0.05555556 0.7222222
4 0.000000 1.61111111 3.5277778 2.11111111 4.5000000
Caecum.trachea Donax.venustus Mysta.picta Glycera.sp. Odostomia.plicata
1 0.000000 0.000000 0.1045752 1.663399 0.2222222
2 1.777778 0.000000 0.1944444 0.000000 0.0000000
3 0.000000 1.088235 0.6143791 0.000000 1.6633987
4 2.861111 0.000000 0.5000000 0.500000 0.0000000
Apseudopsis.ostroumovi Cytharella.costulata Kellia.suborbicularis
1 0.60457516 0.36928105 1.882353
2 0.00000000 0.00000000 0.000000
3 0.05555556 0.05555556 0.000000
4 1.02777778 0.00000000 0.000000
Leptosynapta.inhaerens Spionidae Tritia.reticulata Gastrosaccus.sanctus
1 0.00 2.7222222 0.369281 0.05555556
2 0.00 0.1111111 1.777778 1.25000000
3 0.00 0.0000000 0.000000 0.30065359
4 8.75 1.0277778 0.000000 0.00000000
Micromaldane.ornithochaeta Phoronida Holothuroidea Nereis.perivisceralis
1 0.000000 0.73202614 0.000000 0.00
2 0.000000 0.00000000 0.000000 0.00
3 1.075163 0.05555556 0.000000 0.00
4 0.000000 0.00000000 1.861111 2.75
Phyllodoce.sp. Amphitritides.gracilis Gammaridae Amphipoda Harmothoe.imbricata
1 0.10457516 0.000000 0.00000000 0.000000 2.01634
2 0.44444444 0.000000 2.94444444 5.444444 0.00000
3 0.05555556 0.000000 0.05555556 0.000000 0.00000
4 0.77777778 4.111111 0.00000000 0.000000 0.00000
Mytilus.galloprovincialis Pholoe.inornata Pisione.remota Ampithoe.sp.
1 0.05555556 0.000000 0.05555556 0.0000000
2 0.00000000 0.000000 0.52777778 1.0000000
3 0.00000000 0.000000 0.00000000 0.2647059
4 4.00000000 3.944444 0.19444444 0.0000000
Paradoneis.harpagonea Platyhelminthes Platynereis.dumerilii Rissoa.membranacea
1 0.0000000 0.0000000 0.00 0.0000000
2 0.0000000 0.1944444 0.00 0.0000000
3 0.5882353 0.0000000 0.00 0.9411765
4 0.0000000 0.2777778 2.75 0.0000000
Brachynotus.sexdentatus Eulalia.viridis Microdeutopus.gryllotalpa Nereis.pelagica
1 0.3300654 0.05555556 0.2222222 0.000000
2 0.0000000 0.00000000 0.1111111 0.000000
3 0.0000000 0.00000000 0.1045752 0.000000
4 0.0000000 1.02777778 0.0000000 2.777778
Rapana.venosa Sabellaria.taurica Steromphala.divaricata Aricidea.claudiae
1 0.3300654 0.3300654 0.0000000 0.5359477
2 0.0000000 0.0000000 0.0000000 0.0000000
3 0.0000000 0.0000000 0.0000000 0.0000000
4 0.0000000 0.0000000 0.7777778 0.0000000
Monocorophium.insidiosum Nereis.pulsatoria Parthenina.interstincta Protodrilus.sp.
1 0.0000000 0.26470588 0.05555556 0.000000
2 0.0000000 0.00000000 0.00000000 1.777778
3 0.5359477 0.05555556 0.14705882 0.000000
4 0.0000000 0.00000000 0.00000000 0.000000
Sternaspis.scutata Tellina.fabula Cerastoderma.edule Chondrochelia.savignyi
1 0.000000 0.0000000 0.0000000 0
2 1.777778 0.0000000 0.0000000 0
3 0.000000 0.4183007 0.2647059 0
4 0.000000 0.0000000 0.0000000 1
Crassicorophium.crassicorne Magelona.rosea Polititapes.aureus Polychaeta.larvae
1 0 0.0000000 0.1470588 0.0
2 1 0.0000000 0.0000000 0.5
3 0 0.1470588 0.0000000 0.0
4 0 0.0000000 0.0000000 0.0
Retusa.variabilis Brachystomia.scalaris Eteone.flava Hydrobia.sp. Mactra.stultorum
1 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000
2 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000
3 0.1470588 0.2222222 0.1045752 0.05555556 0.1045752
4 0.0000000 0.0000000 0.0000000 0.11111111 0.0000000
Nephtyidae Papillicardium.papillosum Abra.sp. Acanthocardia.tuberculata Actiniaria
1 0.0000000 0.0000000 0.00000000 0.00000000 0.00000000
2 0.0000000 0.0000000 0.00000000 0.00000000 0.00000000
3 0.2222222 0.1045752 0.05555556 0.05555556 0.05555556
4 0.0000000 0.0000000 0.00000000 0.00000000 0.00000000
Cardiidae Cerastoderma.glaucum Cumopsis.goodsir Elaphognathia.bacescoi
1 0.00000000 0.05555556 0.00000000 0.0000000
2 0.00000000 0.00000000 0.00000000 0.0000000
3 0.05555556 0.00000000 0.05555556 0.0000000
4 0.00000000 0.00000000 0.00000000 0.1111111
Eurydice.pontica Hydrobia.acuta Iphinoe.sp. Lysidice.ninetta Maldanidae
1 0.0000000 0.05555556 0.0000000 0.0000000 0.00000000
2 0.1111111 0.00000000 0.1111111 0.1111111 0.00000000
3 0.0000000 0.00000000 0.0000000 0.0000000 0.05555556
4 0.0000000 0.00000000 0.0000000 0.0000000 0.00000000
Microdeutopus.sp. Neanthes.sp. Nereididae Pestarella.candida Pisces.larvae
1 0.00000000 0.05555556 0.0000000 0.0000000 0.05555556
2 0.00000000 0.00000000 0.1111111 0.0000000 0.00000000
3 0.05555556 0.00000000 0.0000000 0.0000000 0.00000000
4 0.00000000 0.00000000 0.0000000 0.1111111 0.00000000
Retusa.truncatula Rhithropanopeus.harrisii Stenothoe.monoculoides Upogebia.pusilla
1 0.00000000 0.0000000 0.0000000 0.0000000
2 0.00000000 0.0000000 0.0000000 0.0000000
3 0.05555556 0.0000000 0.0000000 0.0000000
4 0.00000000 0.1111111 0.1111111 0.1111111
It’s not perfect, but it’s not too terrible either.
Everything looks more or less fine; fit the model.
glms.lvm.sand <- manyglm(zoo.mvabnd.sand ~ lvm.clusters.sand,
family = "negative.binomial")
Explore the fit (residuals, diagnostic plots, etc.).
png(filename = here(figures.dir, "diag_pl_sand1.png"), width = 16, height = 10, units = "cm", res = 300)
plot.manyglm(glms.lvm.sand, which = 2)
dev.off()
null device
1
I really don’t like the rainbow palette, but I would like to include these plots in my thesis results.. Will have to do something about it, just not right now.
Save the model!
write_rds(glms.lvm.sand,
here(save.dir, "glms_lvm_sand.RDS"))
Let’s see the model summary (NB takes a LOT of time if there are many resamplings!).
(glms.lvm.sand.summary <- summary(glms.lvm.sand,
test = "LR", p.uni = "adjusted",
nBoot = 999, ## limit the number of permutations if you just want to check it out
show.time = "all")
)
The factor (here - groups outlined by the LVM) is highly significant according to the models.
This also allows us to see which species exhibit a response to the chosen factor. The LR (likelihood ratio) statistic is used as a measure of the strength of individual taxon contributions to the observed patterns. I’ll save the summary for safekeeping, but I’ll also run an anova - to get an analysis of deviance table on the model fit (also better for extracting the species contributions, or at least I know how to do it).
write_rds(glms.lvm.sand.summary,
here(save.dir, "glms_lvm_sand_summary.RDS"))
Run the anova on the model.
(glms.lvm.sand.aov <- anova.manyglm(glms.lvm.sand,
test = "LR", p.uni = "adjusted",
nBoot = 999, ## limit the number of permutations for a shorter run time
show.time = "all")
)
I probably shouldn’t have printed all this out, but oh well who cares.
Save the ANOVA, too.
write_rds(glms.lvm.sand.aov,
here(save.dir, "glms_lvm_sand_anova.RDS"))
NOW let’s get the taxa with the highest contributions to the tested pattern (here - clusters in the LVM, which are really the different soft-bottom habitats).
top_n_sp_glm <- function(glms.aov, tot.dev.expl = 0.75) {
## helper retrieving the top n species with the highest contribution to the patterns tested by the GLMs, in decreasing order.
## Arguments: glms.aov - results from an ANOVA on the fitted GLMs
## dev.explained - proportion of explained deviance to use as cutoff
## get the change in deviance due to the tested pattern (= 2nd row from table of univariate test stats), and sort the species in order of decreasing contribution
uni.sorted <- sort(glms.aov$uni.test[2, ], decreasing = TRUE, index.return = FALSE)
## start at 10 species and check how much of the deviance is explained by their contributions. Repeat, increasing by increments of 10 until the desired explained deviance (set at function call) is reached.
top.n.sp <- 10
dev.expl <- sum(uni.sorted[1:top.n.sp])/sum(uni.sorted)
while(dev.expl < tot.dev.expl) {
top.n.sp <- top.n.sp + 10
dev.expl <- sum(uni.sorted[1:top.n.sp])/sum(uni.sorted)
}
## print the total deviance explained - just for information
print(paste("Total deviance explained:", round(dev.expl, 3)))
## return the final top species (and their univariate contributions, just in case)
top.sp <- uni.sorted[1:top.n.sp]
return(top.sp)
}
## get the top contributing species for the initial sand GLMs
(top.sp.glms.lvm.sand <- top_n_sp_glm(glms.lvm.sand.aov, tot.dev.expl = 0.75)
)
## unfortunately, mvabund likes to rename my species when converting the data to matrix (no spaces in names), and since I'm going to look them up in my initial untransformed count data, I have to change them back..
names(top.sp.glms.lvm.sand) <- names(top.sp.glms.lvm.sand) %>%
str_replace(pattern = "\\.", replacement = " ")
top.sp.glms.lvm.sand
Try to plot these top contributing species - for whatever that’s worth, because 50 species on a plot is a monstrosity.
## get the species and their abundances from the original count data, and transform them to long format
(abnd.top.sp.glms.lvm.sand <- zoo.abnd.sand %>%
select(station, names(top.sp.glms.lvm.sand)) %>%
gather(key = "species", value = "count", -station) %>%
## turn species into a factor, or you'll be very very sorry later, when they're out of order on the plot. NB need to be in REVERSE order, because ggplot plots from bottom to top, and I want the top-contributing species on top.
mutate(species = factor(species, levels = rev(names(top.sp.glms.lvm.sand))))
)
plot_top_n <- function(top.n.sp.data, mapping, labs.legend, lab.y, palette) {
## helper for plotting top n species. Was hoping to avoid repeating it from way back when, but no dice.
## Arguments: top.n.sp.data - data frame (long) of top species' counts/biomasses at the different stations
## mapping - mappings of the aesthetics
## labs.legend - labels the use for the legend entries
## lab.y - custom label for y axis
## palette - custom colour palette (for consistency with other plots)
ggplot(top.n.sp.data, mapping) +
geom_point(alpha = 0.75) + # make points larger & partially transparent
scale_color_brewer(palette = palette, labels = labs.legend) +
ylab(lab.y) +
coord_flip()
}
(plot.top.sp.glms.lvm.sand <- plot_top_n(abnd.top.sp.glms.lvm.sand,
mapping = aes(x = species, y = log_y_min(count), colour = station),
labs.legend = paste0("S", as.numeric(unique(abnd.top.sp.glms.lvm.sand$station))),
lab.y = "Abundance (log(y/min + 1))",
palette = "Set2"
) +
theme(legend.position = "top")
)
Well this is a nightmarish plot.. I’ll probably just put this awfulness in a table and call it a day, or play with lvsplot and the modeled ordination plot, if a plot is what’s needed.
Extract the top-contributing species to each cluster (this same nightmare above, but as a table). This chunk is hopelessly ugly and clumsy (and I’ll have to repeat it for the seagrass, too!), but I’m tired of being stuck on this. I still have many, MANY more things to do, and more time-consuming ones too..
top_sp_glms_table <- function(manyglms.obj.smry, group, p = 0.05) {
### extracts the top species in a group for which there is an observed effect in a manyglm test, at the specified probability level.
### Returns: tibble with the top species for the specified group/cluster, sorted (descending) by univariate LR value of the species, significant at the given p level.
## extract the univariate LR coefficients of the species and their p-values
sp_univar <- as_tibble(manyglms.obj.smry$uni.test, rownames = "species")
sp_p <- as_tibble(manyglms.obj.smry$uni.p, rownames = "species")
## combine in the same tibble
sp_all <- left_join(sp_univar, sp_p, by = "species")
## rename the columns
sp_all <- sp_all %>%
rename_at(vars(contains(".x")), list(~str_replace_all(., pattern = ".x", ".LR"))) %>%
rename_at(vars(contains(".y")), list(~str_replace_all(., pattern = ".y", ".p")))
## filter only the group/cluster we want, at the p-level we want
sp_all_flt <- sp_all %>%
select(species, contains(group)) %>%
filter_at(vars(contains(".p")), all_vars(. < p)) %>%
arrange_at(vars(contains(".LR")), list(~desc(.)))
}
top.sp.abnd.glms.lvm.sand <- lapply(names(glms.lvm.sand.summary$aliased), function(x) top_sp_glms_table(glms.lvm.sand.summary, x, p = 0.05))
## fix species names (remove dot)
top.sp.abnd.glms.lvm.sand <- lapply(top.sp.abnd.glms.lvm.sand, function(x) x %>% mutate(species = str_replace(species, pattern = "\\.", replacement = " ")))
## rename columns (= group names) - right now they are something like "lvm.clusters.sand2" etc.
top.sp.abnd.glms.lvm.sand <- lapply(top.sp.abnd.glms.lvm.sand, function(x) x %>% rename_at(vars(contains("lvm.clusters.sand")), list(~str_replace_all(., pattern = "lvm.clusters.sand", "group_"))))
top.sp.abnd.glms.lvm.sand <- lapply(top.sp.abnd.glms.lvm.sand, function(x) x %>% rename_at(vars(contains("Intercept")), list(~str_replace_all(., pattern = "\\(Intercept\\)", "group_1"))))
## pull the abundances from the original count df and add to the summary glm tables
## make a long df of abundances & add clusters
zoo.abnd.sand.long <- zoo.abnd.sand %>%
select(-c(month:replicate)) %>%
gather(key = "species", value = "count", -station) %>%
mutate(group = case_when(station %in% c("Kraimorie", "Chukalya") ~ 1,
station == "Akin" ~ 2,
station %in% c("Sozopol", "Paraskeva") ~ 3,
station == "Agalina" ~ 4))
## sum sp abundances by group; nest by group
zoo.abnd.sand.long.smry <- zoo.abnd.sand.long %>%
group_by(species, group) %>%
summarise(total_count = sum(count)) %>%
group_by(group) %>%
nest()
## add the counts to the group dfs - wow that's an ugly, ugly hack. Wish I had more time to write this up properly..
top.sp.abnd.glms.lvm.sand <- map2(top.sp.abnd.glms.lvm.sand, zoo.abnd.sand.long.smry %>% pull(group), ~left_join(.x, zoo.abnd.sand.long.smry %>% filter(group == .y) %>% unnest(), by = "species"))
## since these are sum counts over all the replicates (that's why the monstrous numbers), average them to be mean counts per group. NB different groups consist of different numbers of replicates, b.c. some groups consist of more than one station
(top.sp.abnd.glms.lvm.sand <- map2(top.sp.abnd.glms.lvm.sand, c(18, 9, 18, 9), function(x, y) x %>% mutate(mean_count = total_count/y))
)
To determine the relative taxon contribution to patterns: LR statistic - a measure of strength of individual taxon contributions. LR expresses how many times more likely the data are under one model than the other. This likelihood ratio, or equivalently its logarithm, can then be used to compute a p-value, or, compared to a critical value, to decide whether to reject the null model in favour of the alternative model.
In this case, the model shows which species exhibit a reaction based on the chosen groups - in other words, which species are more likely to be more/less abundant in each group.
For group 1 (= S1-S2), the species/taxa with significantly higher abundance are: Oligochaeta, H. filiformis, P. kefersteini, M. palmata, P. cirrifera, A. diadema (among others); and the ones with significantly lower abundance - even 0, in some cases - S. bidentata, B.lanceolatum, M. papillicornis, Melita palmata, P. jubatus, and so on.
For group 2 (= S3), the species with higher abundance are: B. lanceolatum, O. limacina, Oligochaeta (this is this strange artifact of 2013), P. kefersteini, L. flavocapitatus. The species with lower abundance are: H. filiformis, A. kagoshimensis, M. stammeri, Melinna palmata, etc. For group 3 (= S4-S6), the species with higher abundance are: C. gallina, L. mediterraneum - with very high dominance over practically all others; also Pseudocuma longicorne, Spio filicornis. The species with lower abundance are: H. filiformis, Oligochaetes (to a certain extent - they are still present, though), A. kagoshimensis, L. koreni, Harmothoe reticulata, Iphinoe tenella, Leiochone leiopygos.
For group 4 (= S5), the species with higher abundance are: Microdeutopus versiculatus, Eurydice dollfusi, Melita palmata, Polygordius neapolitanus, Polycirrus caliendrum, Polycirrus jubatus, Streptosyllis bidentata. The species with lower abundance are: A. kagoshimensis, Melinna palmata, P. cirrifera, P. ciliata, A. alba, I. tenella.
I love how the species with the highest variances (e.g. C. gallina, the most conspicuous example) are consistently pushed back - have lower LR scores. This is very good - C. gallina in particular is dominant in group 3, but is present also in all other groups - its substrate/depth preferences are very wide, so this is not uncommon. It’s not automatically pushed to the top of the list, but its reaction is detected by the manyGLM test. Neat! Contrast to the SIMPER results, where the species with the highest variance are consistently at the top - they contribute the most to the similarity, as per the test definition.
I’m going to save these as separate files (manually), then format them as tables - I know it’s a shame, but I’m too frustrated to figure out how to do it programmatically.
I’ll also put them in a word table in my final text, because I don’t want to deal with a million separate ones (embedded excel tables don’t split over multiple pages).
NB In my text, I’m switching the names/places of group 3 and 4, to be consistent with the SIMPER groups (I’m NOT going to repeat all this just to have the numbers match up). So the file names, table names, etc. remain as above. But in the text, I’ll have the following: group 1 = S1-S2, group 2 = S3, group 3 = S5, group 4 = S4-S6. REMEMBER THIS SO THERE IS NO CONFUSION!
Now, let’s try to see a different thing - which environmental parameters best describe the species response.
I’m going to use the PCA-filtered environmental data - it’s still going to be a slog, with 7 potential predictors..
First, construct the formula for the model - will do it separately in case I need to update it later, etc. This is the full formula with all explanatory variables.
(formula.env.glms.sand <- formula(paste("zoo.mvabnd.sand ~",
paste(env.sand %>% select(-station) %>% names(), collapse = "+")))
)
Fit the GLMs to the sand abundance data.
env.glms.sand <- manyglm(formula.env.glms.sand,
data = env.sand,
family = "negative.binomial")
Explore the fit (residuals, diagnostic plots, etc.).
## residuals vs fitted values
plot(env.glms.sand)
## all traditional (g)lm diagnostic plots
plot.manyglm(env.glms.sand, which = 1:3)
# ### source mvabund GLM plotting functions modified to use a grey palette - I just can't redo these plots on my own, the function is doing too complicated things internally to scale the x and y axes
# source(here(functions.dir, "default.plot.manyglm_grey.R"))
# source(here(functions.dir, "plot.manyglm_grey.R"))
#
# par(mfrow = c(2,2))
# lapply(1:3, function(i) plot.manyglm.grey(glms.lvm.sand, which = i, sub.caption = ""))
# par(mfrow = c(1, 1))
Well, it’s good enough if you ask me (still the kinda strange “line” at lin.pred = -6; otherwise residuals are random enough).
Save the model!
write_rds(env.glms.sand,
here(save.dir, "glms_env_sand.RDS"))
Before anything else, I want to try and reduce the model a little - to improve the fit/reduce run time.
The automatic step functions that eliminate/add model terms sequentially don’t work - they fail at the last step with a cryptic error about differing numbers of rows - I assume because manyglm has as left side term the whole community abundance matrix, and the functions don’t really know how to deal with that. I don’t understand enough about their internals to fix the problem, so I’m just going to write my own little automation based on the function drop1.
evaluate_glms_env <- function(full.mod) {
### sequentially eliminate model terms in manyglms vs environmental parameters, and find the best model based on lowest AIC score.
### Arguments: full.mod - full model fit
### Returns: best manyglm model of environmental parameters; prints out the best model formula
### Dependencies: tidyverse
## get the starting formula (= full model with all variables)
start.formula <- formula(full.mod)
drop_var <- function(mod) {
### helper picking the next variable to drop from a model to improve the fit (based on AIC)
## check the model AICs if variables are dropped one by one
drop1.df <- as_tibble(drop1(mod), rownames = "drop_var") %>% arrange(AIC)
## pick the variable to drop next - the one resulting in the largest decrease in AIC
drop.var <- drop1.df %>% filter(AIC == min(AIC)) %>% pull(drop_var)
return(drop.var)
}
## pick the variable to drop next
drop.var <- drop_var(full.mod)
if(drop.var != "<none>") {
## update the model formula, dropping the variable resulting in the largest decrease in AIC; then apply it to the model.
new.formula <- update.formula(start.formula, paste0("~. -", drop.var))
new.mod <- update(full.mod, new.formula)
## identify a new variable to drop that lowers the AIC
drop.var <- drop_var(new.mod)
## repeat the steps above until the function can no longer find such a variable (i.e., dropping more variables doesn't improve the model fit)
while(drop.var != "<none>") {
new.formula <- update.formula(new.formula, paste0("~. -", drop.var))
new.mod <- update(full.mod, new.formula)
drop.var <- drop_var(new.mod)
}
## print out the best model formula
print(paste("Best model: ", paste(deparse(new.formula), collapse = "")))
return(new.mod)
} else {
## if the starting model is the best, print its formula (fat chance!)
print(paste("Best model: ", paste(deparse(start.formula), collapse = "")))
return(full.mod)
}
}
Select the best reduced model of environmental variables for the sand stations.
## selection function defined in the sand section
top.env.glm.red.sand <- evaluate_glms_env(env.glms.sand)
Check its fit.
## residuals vs fitted values
plot(top.env.glm.red.sand)
## all traditional (g)lm diagnostic plots
plot.manyglm(top.env.glm.red.sand, which = 1:3)
I think it’s fine; might even be better than the full model.. Save it, too.
write_rds(top.env.glm.red.sand,
here(save.dir, "glms_top_env_red_sand.RDS"))
Run ANOVA on this model.
(top.env.glm.red.sand.aov <- anova.manyglm(top.env.glm.red.sand,
test = "LR", p.uni = "adjusted",
nBoot = 999, ## limit the number of permutations for a shorter run time
show.time = "all")
)
Resampling begins for test 1.
Resampling run 0 finished. Time elapsed: 0.01 minutes...
Resampling run 100 finished. Time elapsed: 0.58 minutes...
Resampling run 200 finished. Time elapsed: 1.18 minutes...
Resampling run 300 finished. Time elapsed: 1.76 minutes...
Resampling run 400 finished. Time elapsed: 2.35 minutes...
Resampling run 500 finished. Time elapsed: 2.95 minutes...
Resampling run 600 finished. Time elapsed: 3.56 minutes...
Resampling run 700 finished. Time elapsed: 4.18 minutes...
Resampling run 800 finished. Time elapsed: 4.79 minutes...
Resampling run 900 finished. Time elapsed: 5.43 minutes...
Resampling begins for test 2.
Resampling run 0 finished. Time elapsed: 0.00 minutes...
Resampling run 100 finished. Time elapsed: 0.36 minutes...
Resampling run 200 finished. Time elapsed: 0.73 minutes...
Resampling run 300 finished. Time elapsed: 1.09 minutes...
Resampling run 400 finished. Time elapsed: 1.47 minutes...
Resampling run 500 finished. Time elapsed: 1.83 minutes...
Resampling run 600 finished. Time elapsed: 2.21 minutes...
Resampling run 700 finished. Time elapsed: 2.58 minutes...
Resampling run 800 finished. Time elapsed: 2.95 minutes...
Resampling run 900 finished. Time elapsed: 3.30 minutes...
Resampling begins for test 3.
Resampling run 0 finished. Time elapsed: 0.00 minutes...
Resampling run 100 finished. Time elapsed: 0.35 minutes...
Resampling run 200 finished. Time elapsed: 0.70 minutes...
Resampling run 300 finished. Time elapsed: 1.04 minutes...
Resampling run 400 finished. Time elapsed: 1.38 minutes...
Resampling run 500 finished. Time elapsed: 1.73 minutes...
Resampling run 600 finished. Time elapsed: 2.08 minutes...
Resampling run 700 finished. Time elapsed: 2.43 minutes...
Resampling run 800 finished. Time elapsed: 2.79 minutes...
Resampling run 900 finished. Time elapsed: 3.15 minutes...
Resampling begins for test 4.
Resampling run 0 finished. Time elapsed: 0.00 minutes...
Resampling run 100 finished. Time elapsed: 0.32 minutes...
Resampling run 200 finished. Time elapsed: 0.64 minutes...
Resampling run 300 finished. Time elapsed: 0.95 minutes...
Resampling run 400 finished. Time elapsed: 1.27 minutes...
Resampling run 500 finished. Time elapsed: 1.58 minutes...
Resampling run 600 finished. Time elapsed: 1.89 minutes...
Resampling run 700 finished. Time elapsed: 2.20 minutes...
Resampling run 800 finished. Time elapsed: 2.51 minutes...
Resampling run 900 finished. Time elapsed: 2.83 minutes...
Time elapsed: 0 hr 16 min 18 sec
Analysis of Deviance Table
Model: manyglm(formula = zoo.mvabnd.sand ~ PO4 + seston + LUSI + gravel,
Model: family = "negative.binomial", data = env.sand)
Multivariate test:
Res.Df Df.diff Dev Pr(>Dev)
(Intercept) 53
PO4 52 1 1147.2 0.001 ***
seston 51 1 816.4 0.001 ***
LUSI 50 1 728.5 0.001 ***
gravel 49 1 1128.7 0.001 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Univariate Tests:
Abra.alba Abra.sp. Acanthocardia.tuberculata
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
PO4 1.011 1.000 0.513 1.000 3.576 0.979
Acari Actiniaria Alitta.succinea
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
PO4 10.498 0.120 0.513 1.000 0.065 1.000
Ampelisca.diadema Amphibalanus.improvisus Amphipoda
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
PO4 12.592 0.032 11.462 0.081 0.552
Amphitritides.gracilis Ampithoe.sp.
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
PO4 1.000 2.571 1.000 0.032 1.000
Anadara.kagoshimensis Aonides.paucibranchiata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
PO4 67.791 0.001 5.965 0.659
Apseudopsis.ostroumovi Aricidea.claudiae
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
PO4 0.375 1.000 7.633 0.402
Bathyporeia.guilliamsoniana Bittium.reticulatum
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
PO4 0.002 1.000 22.784 0.002
Bodotria.arenosa Brachynotus.sexdentatus
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
PO4 1.3 1.000 11.298 0.085
Brachystomia.scalaris Branchiostoma.lanceolatum
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
PO4 0.696 1.000 0.057 1.000
Caecum.armoricum Caecum.trachea Capitella.capitata
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
PO4 7.463 0.427 1.119 1.000 1.193
Capitella.minima Cardiidae Cerastoderma.edule
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
PO4 1.000 2.53 1.000 0.513 1.000 7.708
Cerastoderma.glaucum Chamelea.gallina
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
PO4 0.384 1.415 1.000 27.316 0.001
Chondrochelia.savignyi Crassicorophium.crassicorne
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
PO4 1.246 1.000 0.459 1.000
Cumopsis.goodsir Cytharella.costulata Decapoda.larvae
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
PO4 3.576 0.969 16.556 0.008 4.389
Dinophilus.gyrociliatus Diogenes.pugilator
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
PO4 0.845 3.462 0.988 2.328 1.000
Donax.venustus Elaphognathia.bacescoi Eteone.flava
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
PO4 11.663 0.069 1.021 1.000 7.158
Eteone.sp. Eulalia.viridis Eunice.vittata
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
PO4 0.480 3.247 0.996 0.508 1.000 5.957
Eurydice.dollfusi Eurydice.pontica
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
PO4 0.659 7.553 0.419 0.243 1.000
Exogone.naidina Gammaridae Gastrosaccus.sanctus
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
PO4 1.544 1.000 0.215 1.000 0.001 1.000
Genetyllis.tuberculata Glycera.sp. Glycera.tridactyla
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
PO4 9.748 0.149 1.819 1.000 4.038
Harmothoe.imbricata Harmothoe.reticulata
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
PO4 0.882 7.572 0.416 6.871 0.582
Heteromastus.filiformis Hirudinea Holothuroidea
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
PO4 55.658 0.001 4.073 0.876 6.327 0.634
Hydrobia.acuta Hydrobia.sp. Iphinoe.sp.
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
PO4 3.566 0.988 3.379 0.993 0.243 1.000
Iphinoe.tenella Kellia.suborbicularis Lagis.koreni
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
PO4 38.241 0.001 12.137 0.052 48.067
Leiochone.leiopygos Lentidium.mediterraneum
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
PO4 0.001 20.554 0.003 16.232 0.009
Leptosynapta.inhaerens Lindrilus.flavocapitatus
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
PO4 3.875 0.905 4.746 0.795
Loripes.orbiculatus Lucinella.divaricata Lysidice.ninetta
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
PO4 14.871 0.012 19.259 0.004 0.243
Mactra.stultorum Magelona.papillicornis
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
PO4 1.000 2.401 1.000 19.059 0.004
Magelona.rosea Maldanidae Melinna.palmata
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
PO4 4.77 0.795 3.576 0.979 75.355 0.001
Melita.palmata Microdeutopus.gryllotalpa Microdeutopus.sp.
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
PO4 14.255 0.019 0.001 1.000 3.576
Microdeutopus.versiculatus Micromaldane.ornithochaeta
Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
PO4 0.971 9.28 0.172 10.809
Micronephthys.stammeri Microphthalmus.fragilis
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
PO4 0.098 0.016 1.000 0.028 1.000
Microphthalmus.sp. Monocorophium.acherusicum
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
PO4 8.692 0.220 3.627 0.966
Monocorophium.insidiosum Mysta.picta Mytilaster.lineatus
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
PO4 7.646 0.401 4.088 0.876 9.876
Mytilus.galloprovincialis Neanthes.sp. Nemertea
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
PO4 0.143 0.507 1.000 1.415 1.000 0.718
Nephtyidae Nephtys.cirrosa Nereididae
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
PO4 1.000 3.733 0.931 3.455 0.989 0.243 1.000
Nereis.pelagica Nereis.perivisceralis Nereis.pulsatoria
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
PO4 1.278 1.000 5.021 0.764 1.04
Nototropis.guttatus Odostomia.plicata Oligochaeta
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
PO4 1.000 0.017 1.000 2.158 1.000 12.804
Ophelia.limacina Papillicardium.papillosum
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA> <NA>
PO4 0.032 1.13 1.000 1.026 1.000
Paradoneis.harpagonea Parthenina.interstincta
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA>
PO4 16.021 0.009 0.628 1.000
Parvicardium.exiguum Perinereis.cultrifera
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA>
PO4 14.43 0.013 6.906 0.582
Perioculodes.longimanus Pestarella.candida Pholoe.inornata
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept) <NA> <NA> <NA> <NA> <NA>
PO4 13.283 0.026 1.021 1.000 2.532
Phoronida Phyllodoce.sp. Pisces.larvae
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA> <NA> <NA> <NA>
PO4 1.000 14.128 0.020 0.214 1.000 3.566 0.979
Pisione.remota Pitar.rudis Platyhelminthes
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA> <NA> <NA>
PO4 0.184 1.000 2.33 1.000 1.258 1.000
Platynereis.dumerilii Polititapes.aureus Polychaeta.larvae
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept) <NA> <NA> <NA> <NA> <NA>
PO4 2.51 1.000 5.232 0.758 0.667
Polycirrus.caliendrum Polycirrus.jubatus
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA> <NA>
PO4 1.000 8.778 0.213 14.046 0.021
Polydora.ciliata Polygordius.neapolitanus
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA>
PO4 33.47 0.001 1.738 1.000
Prionospio.cirrifera Protodorvillea.kefersteini
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA>
PO4 48.371 0.001 2.268 1.000
Protodrilus.sp. Pseudocuma.longicorne Rapana.venosa
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept) <NA> <NA> <NA> <NA> <NA>
PO4 0.497 1.000 13.533 0.024 16.327
Retusa.truncatula Retusa.variabilis
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA> <NA>
PO4 0.008 0.513 1.000 4.77 0.795
Rhithropanopeus.harrisii Rissoa.membranacea
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA>
PO4 1.021 1.000 2.111 1.000
Sabellaria.taurica Salvatoria.clavata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA>
PO4 16.327 0.008 2.946 0.999
Schistomeringos.rudolphi Sphaerosyllis.hystrix
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA>
PO4 11.176 0.087 3.042 0.998
Spio.filicornis Spionidae Spisula.subtruncata
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA> <NA> <NA>
PO4 8.141 0.262 0.278 1.000 34.272 0.001
Stenothoe.monoculoides Sternaspis.scutata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA>
PO4 1.021 1.000 0.497 1.000
Steromphala.divaricata Streptosyllis.bidentata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA>
PO4 3.62 0.966 14.133 0.020
Syllis.gracilis Syllis.hyalina Tellina.fabula
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA> <NA> <NA>
PO4 1.059 1.000 11.199 0.086 7.846 0.291
Tellina.tenuis Tritia.neritea Tritia.reticulata
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA> <NA> <NA>
PO4 5.219 0.758 9.738 0.149 7.578 0.416
Turbellaria Upogebia.pusilla
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA>
PO4 0.973 1.000 1.021 1.000
[ reached getOption("max.print") -- omitted 3 rows ]
Arguments:
Test statistics calculated assuming uncorrelated response (for faster computation)
P-value calculated using 999 resampling iterations via PIT-trap resampling (to account for correlation in testing.
So, it turns out that the long-term water column eutrophication (PO4, seston), the anthropogenic pressure in general (LUSI), and the sediment composition (gravel) explain the observed community patterns best.
Save the ANOVA - I really, really don’t want to have to repeat it.
write_rds(top.env.glm.red.sand.aov,
here(save.dir, "glms_top_env_red_sand_anova.RDS"))
Get the taxa with the highest contributions to the tested pattern (here - species most affected by changes in water/environmental quality parameters).
## get the top contributing species for the environmental parameter sand GLMs
(top.sp.glms.env.red.sand <- top_n_sp_glm(top.env.glm.red.sand.aov, tot.dev.expl = 0.75)
)
[1] "Total deviance explained: 0.817"
Melinna.palmata Anadara.kagoshimensis Heteromastus.filiformis
75.354817 67.790515 55.657621
Prionospio.cirrifera Lagis.koreni Iphinoe.tenella
48.371413 48.067236 38.240757
Spisula.subtruncata Polydora.ciliata Chamelea.gallina
34.272034 33.469990 27.316472
Bittium.reticulatum Leiochone.leiopygos Lucinella.divaricata
22.783736 20.553990 19.259453
Magelona.papillicornis Cytharella.costulata Sabellaria.taurica
19.059044 16.556175 16.326603
Rapana.venosa Lentidium.mediterraneum Paradoneis.harpagonea
16.326603 16.231882 16.021000
Loripes.orbiculatus Parvicardium.exiguum Melita.palmata
14.870959 14.429942 14.255372
Streptosyllis.bidentata Phoronida Polycirrus.jubatus
14.133303 14.128398 14.045827
Pseudocuma.longicorne Perioculodes.longimanus Oligochaeta
13.533423 13.282665 12.804148
Ampelisca.diadema Kellia.suborbicularis Donax.venustus
12.591581 12.137452 11.662509
Amphibalanus.improvisus Brachynotus.sexdentatus Syllis.hyalina
11.462125 11.298400 11.198693
Schistomeringos.rudolphi Micromaldane.ornithochaeta Acari
11.176121 10.809322 10.498478
Mytilaster.lineatus Genetyllis.tuberculata Tritia.neritea
9.875996 9.747742 9.737885
Microdeutopus.versiculatus Polycirrus.caliendrum Microphthalmus.sp.
9.280497 8.778395 8.692123
Spio.filicornis Tellina.fabula Cerastoderma.edule
8.140655 7.845543 7.707880
Monocorophium.insidiosum Aricidea.claudiae Tritia.reticulata
7.646026 7.632964 7.577649
Harmothoe.imbricata Eurydice.dollfusi
7.571700 7.553233
## unfortunately, mvabund likes to rename my species when converting the data to matrix (no spaces in names), and since I'm going to look them up in my initial untransformed count data, I have to change them back.. DON'T BE IN A HURRY TO DO THAT IF YOU WANT TO SUBSET THE ORIGINAL MATRIX BEFORE RUNNING TRAITGLM
names(top.sp.glms.env.red.sand) <- names(top.sp.glms.env.red.sand) %>%
str_replace(pattern = "\\.", replacement = " ")
I’m going to plot these top contributing species, but I’m not using the plot. At least this time it’s more manageable, but still not presentable enough..
## get the species and their abundances from the original count data, and transform them to long format
abnd.top.sp.glms.env.red.sand <- zoo.abnd.sand %>%
select(station, names(top.sp.glms.env.red.sand)) %>%
gather(key = "species", value = "count", -station) %>%
## turn species into a factor, or you'll be very very sorry later, when they're out of order on the plot. NB need to be in REVERSE order, because ggplot plots from bottom to top, and I want the top-contributing species on top.
mutate(species = factor(species, levels = rev(names(top.sp.glms.env.red.sand)))) %>%
## add clusters from LVM as a column
mutate(group = case_when(station %in% c("Kraimorie", "Chukalya") ~ 1,
station == "Akin" ~ 2,
station %in% c("Sozopol", "Paraskeva") ~ 3,
station == "Agalina" ~ 4))
## add the significant environmental parameters from the model
(abnd.top.sp.glms.env.red.sand <- left_join(abnd.top.sp.glms.env.red.sand,
env.sand %>% select(station, PO4, seston, LUSI, gravel),
by = "station")
)
(plot.top.sp.glms.env.red.sand <- plot_top_n(abnd.top.sp.glms.env.red.sand,
mapping = aes(x = species, y = log_y_min(count), colour = factor(group)),
labs.legend = unique(abnd.top.sp.glms.env.red.sand$group),
lab.y = "Abundance (log(y/min + 1))",
palette = "Set2") +
theme(legend.position = "top")
)
Try another type of chart - to start, a parallel coordinates chart, maybe.
# You need the MASS library
# Vector color
my_colors=colors()[as.numeric(unique(abnd.top.sp.glms.env.red.sand$species))]
# Make the graph !
MASS::parcoord(abnd.top.sp.glms.env.red.sand %>% select(PO4:gravel), col = my_colors)
This is terrible, not to mention non-informative..
Extract the taxon information (univariate tests) from the model ANOVA to present as a table (probably better than this plot, although it’s informative).
## extract the univariate test coefficients (LR) from the environmental model ANOVA. NB keep the row names when converting the matrix to tibble!
table.top.sp.glms.env.red.sand <- as_tibble(top.env.glm.red.sand.aov$uni.test, rownames = "var")
## fix the species names - remove first dor
names(table.top.sp.glms.env.red.sand) <- names(table.top.sp.glms.env.red.sand) %>%
str_replace(pattern = "\\.", replacement = " ")
## subset only the top species (explaining ~75% of the dataset variation)
table.top.sp.glms.env.red.sand <- table.top.sp.glms.env.red.sand %>%
select(var, names(top.sp.glms.env.red.sand))
## transpose, because a table with 50 columns is just unreadable
(table.top.sp.glms.env.red.sand <- table.top.sp.glms.env.red.sand %>%
gather(key = species, value = value, -var) %>%
spread(key = var, value = value) %>%
## arrange as before (terms in the order they appear in the model, and by descending value of the LR for the first model term - here, PO4). Also get rid of the intercept (it's all-NA anyway).
select(species, PO4, seston, LUSI, gravel) %>%
arrange(desc(PO4))
)
Save this to a file - will have to format it as a nice table by hand, unfortunately.
write_csv(table.top.sp.glms.env.red.sand,
here(save.dir, "taxa_contrib_glms_top_env_red_sand.csv"))
Final analysis to try: which species respond differently to different environmental parameters? (= traits analysis - fit single predictive model for all species at all sites, but w/o attempting to explain the different responses using traits - the species ID is used in place of a traits matrix).
NB only use the top species that exhibited a reaction in the environmental model fit (= the ones accounting for ~75% of the total variability), and only the significant predictors - to improve run times.
sp.response.glms.env.red.sand <- traitglm(L = mvabund(zoo.abnd.flt.sand[, names(top.sp.glms.env.red.sand)]),
R = as.matrix(env.sand %>% select(PO4, seston, LUSI, gravel)),
method = "manyglm")
No traits matrix entered, so will fit SDMs with different env response for each spp
sp.response.glms.env.red.sand$fourth.corner
PO4 seston LUSI gravel
names.L.Ampelisca.diadema 3.9255509 -2.9690187 -0.817892764 -0.54638276
names.L.Amphibalanus.improvisus 4.0226904 -3.2428168 -0.613822194 -0.55947227
names.L.Anadara.kagoshimensis 11.2652127 -7.5624621 -2.881581827 -0.56229429
names.L.Aricidea.claudiae 6.4545107 -6.0030102 0.075837116 -0.81658105
names.L.Bittium.reticulatum 3.6728293 -3.2566010 -0.227682651 -0.60908662
names.L.Brachynotus.sexdentatus 12.2714488 -8.7953746 -2.768356521 -1.15909023
names.L.Cerastoderma.edule 86.3144545 -68.5410436 -16.124587996 -15.76009706
names.L.Chamelea.gallina 4.9735567 -4.0806520 -0.915304631 -0.90870404
names.L.Cytharella.costulata 5.9066074 -4.5531016 -1.156170995 -0.64061054
names.L.Donax.venustus 29.2665146 -25.3133054 -3.196472244 -7.48930781
names.L.Eurydice.dollfusi -145.3134807 78.0140182 56.500452934 -1.12048695
names.L.Genetyllis.tuberculata 4.9508698 -4.4226734 -0.717237975 -0.50545895
names.L.Harmothoe.imbricata 8.7688494 -8.0381252 -0.174424875 -1.40001396
names.L.Heteromastus.filiformis 16.4489612 -11.8603189 -4.679979720 -0.53453476
names.L.Iphinoe.tenella 9.6429106 -6.3711957 -2.467970830 -0.34728546
names.L.Kellia.suborbicularis 18.1479826 -13.1444944 -4.028384566 -2.38364695
names.L.Lagis.koreni 5.2752482 -3.9409882 -1.063633817 -0.49005707
names.L.Leiochone.leiopygos 9.4458023 -7.0107877 -2.426883847 -0.49499547
names.L.Lentidium.mediterraneum 2.2184223 -7.7274293 5.119259705 -3.80721597
names.L.Loripes.orbiculatus 11.9138363 -10.5152606 -1.772124177 -1.53849109
names.L.Lucinella.divaricata 4.3089522 -3.8383729 -0.941760570 -0.90040958
names.L.Magelona.papillicornis -0.8618886 -1.6184538 2.126924166 -1.13824665
names.L.Melinna.palmata 87.5950547 -62.3176332 -28.320862498 -0.53785792
names.L.Melita.palmata -1.4068078 -1.4310911 2.175555342 0.04534594
names.L.Microdeutopus.versiculatus 4.5619062 -3.3188059 -1.480537959 0.01704370
names.L.Micromaldane.ornithochaeta 2.0293951 -3.6054716 1.146677297 -1.10488159
names.L.Microphthalmus.sp. 2.4219932 -2.3771910 -0.331905804 -0.45014535
names.L.Monocorophium.insidiosum 91.2704372 -72.4819591 -16.997678796 -16.75142841
names.L.Mytilaster.lineatus 3.6452068 -3.0310816 -0.453771212 -0.56603117
names.L.Oligochaeta 1.1142783 -0.9168759 0.182815777 -0.26804214
names.L.Paradoneis.harpagonea 44.6268239 -34.6189865 -10.110156651 -6.18181742
names.L.Parvicardium.exiguum 5.4679515 -4.2901258 -1.071674867 -0.60043740
names.L.Perioculodes.longimanus 3.8233980 -3.1155895 -0.537523634 -0.55747087
names.L.Phoronida 4.3674755 -4.1186121 -0.107432136 -0.69239101
names.L.Polycirrus.caliendrum -14.6362720 -14.5110317 24.537891679 -3.49252794
names.L.Polycirrus.jubatus -2.4679942 -2.7467481 4.271575754 -0.36118532
names.L.Polydora.ciliata 5.8538836 -4.5251508 -1.054314612 -0.53519651
names.L.Prionospio.cirrifera 4.6312252 -3.5794420 -0.682191827 -0.61411757
names.L.Pseudocuma.longicorne -184.7247688 94.1777708 77.116820970 -2.87857061
names.L.Rapana.venosa 6.2536074 -5.8242086 0.109640812 -0.74952507
names.L.Sabellaria.taurica 6.4201209 -5.9779189 0.094454203 -0.79630833
names.L.Schistomeringos.rudolphi 8.9818953 -6.6900990 -3.539271831 -0.05466283
names.L.Spio.filicornis 23.9025062 -18.4529248 -3.775190226 -6.60866099
names.L.Spisula.subtruncata 8.3692147 -6.4982724 -1.857417192 -1.49042998
names.L.Streptosyllis.bidentata -1.7193139 -2.0772085 3.037260251 -0.16308400
names.L.Syllis.hyalina -1.4595610 -2.6762110 3.349712166 -0.32574204
names.L.Tellina.fabula 91.2704372 -72.4819591 -16.997678816 -16.75142841
names.L.Tritia.neritea -3.8087181 -2.4272915 5.454585358 -1.65134985
names.L.Tritia.reticulata 2.6224780 -1.4687956 0.003212383 -0.74189472
# plot this
a <- max(abs(sp.response.glms.env.red.sand$fourth.corner))
colort <- colorRampPalette(c("blue","white","red"))
plot.spp <- lattice::levelplot(t(as.matrix(sp.response.glms.env.red.sand$fourth.corner)), xlab = "Environmental Variables",
ylab = "Species", col.regions = colort(100), at = seq(-a, a, length = 100),
scales = list(x = list(rot = 45)))
print(plot.spp)
When using LASSO (method = “glm1path”), the algorithm fails to converge - I’m not sure how to interpret it.. Maybe because the function tests each individual species:env.parameter interaction (does it really??), and none of them by themselves are sufficient to explain a species’ response. Not to mention the fact that the samples are not really independent (they are replicates at 6 sites, repeated 3 times).
When using method = “manyglm”, the result is the one shown above. It’s still a bitch to interpret - for example, what is the interpretation of an increase in abundance with for ex. high PO4, but low LUSI? Where are these conditions ever met?
In fact, everything points towards the conclusion that a species response is determined by a combination of eutrophication parameters in its environment (water column characteristics), and the composition of the sediments (organic matter and granulometry).
This is actually sort of similar to the PERMANOVA results, in this particular case. However, it’s much more parsimonious.
In the future, I’m leaning more towards the modeling approach - it allows you to check the model fit to one’s real data; also, there are no data reductions due to calculation of distance matrices.
Import zoobenthic abundance data (cleaned and prepared).
zoo.abnd.zostera <- read_csv(here(save.dir, "abnd_zostera_orig_clean.csv"))
## convert station to factor (better safe than sorry later, when the stations are not plotted in the order I want them)
(zoo.abnd.zostera <- zoo.abnd.zostera %>%
mutate(station = factor(station, levels = c("Poda", "Otmanli", "Vromos", "Gradina", "Ropotamo")))
)
Remove the all-0 species (= not present in the current dataset).
Maybe also remove the singletons (species appearing only once in the whole dataset and represented by a single individual = so rare that it’s unlikely they carry important information, but it would probably improve the run times).
(zoo.abnd.flt.zostera <- zoo.abnd.zostera %>%
select(-c(station:replicate)) %>%
select(which(colSums(.) > 0))
)
Perform a model-based unconstrained ordination by fiting a pure latent variable model (package boral - Hui et al., 2014). This will allow to visualize the multivariate stations x species data - similar to nMDS, can be interpreted in the same way.
I’m including a (fixed) row effect to account for differences in site total abundance - this way, the ordination is in terms of species composition.
NB this takes about a million years to run!
lvm.zostera <- boral(y = zoo.abnd.flt.zostera,
family = "negative.binomial",
## we want to control for site effects - there are 6 sites with 9 replicates each
row.eff = "fixed", row.ids = matrix(rep(1:5, times = c(8, 8, 4, 8, 4)), ncol = 1),
## 2 latent variables = 2 axes on which to represent the zoobenthic data
lv.control = list(num.lv = 2)
# ## example control structure, to check if function does what I want, because otherwise it takes an intolerably long time, and I'll shoot myself if I have to wait for it again
# mcmc.control = list(n.burnin = 10, n.iteration = 100,
# n.thin = 1)
#
)
Check the summary and diagnostic plots for the LVM.
summary(lvm.zostera)
## model fit diagnostic plots
plot(lvm.zostera)
The residuals plots look fine (no patterns in the residuals vs fitted, so variance is homogeneous, the quantile plot shows a (more or less) normal distribution of the residuals) - the model fits the data pretty well.
Save the zostera LVM.
write_rds(lvm.zostera,
here(save.dir, "lvm_zostera.RDS"))
Examine the biplot obtained by fitting the LVM, as well as the 20 most “important” species.
lvsplot(lvm.zostera, jitter = T, biplot = TRUE, ind.spp = 20)
All in all, the final result resembles the nMDS ordination very much - same stretched clusters (Poda + Otmanli, Vromos pretty much apart, Gradina +- Ropotamo). I don’t see much difference with the nMDS. The main difference seems to be the distance between the 2 years for Poda ana Otmanli - the LVM enlarges it. Have to remember to test for year effect! The run time is actually not that bad for the seagrasses. The species singled out as significant are probably somewhat different - have to check!
Redo the biplot, because this one is not very pretty. I’m not adding the species on top, first because I’m too lazy to figure out the procedure for ordering them, and second because the plot gets too busy.
## extract the LV coordinates of the stations from the model, so that the plot can be redone in ggplot
lvs.coord.zostera <- as_tibble(lvm.zostera$lv.median)
## add the stations from the original zoobenthic table (order was not modified)
(lvs.coord.zostera <- lvs.coord.zostera %>%
bind_cols(zoo.abnd.zostera %>% select(station))
)
Make the plot and save it.
(plot.lvm.zostera <- ggplot(lvs.coord.zostera) +
geom_point(aes(x = lv1, y = lv2, colour = station)) +
scale_color_brewer(palette = "Set2", name = "station",
labels = paste0("Z", as.numeric((unique(lvs.coord.zostera %>% pull(station)))))) +
labs(x = "LV1", y = "LV2")
)
Well, this is a weird one - this plot is flipped around 0 compared to the one that boral’s plotting function gives. Otherwise nothing changes - the spatial relationships between samples are preserved. I suppose it doesn’t matter much - the axes are arbitrary after all, but strange that it happens.
## save the LVM plot for the seagrass
ggsave(file = here(figures.dir, "lvm_zostera.png"),
plot.lvm.zostera,
width = 15, units = "cm", dpi = 300)
Fit GLMs to the sites x species matrix to try and explain the observed differences in community structure by the variation of the environmental parameters.
These functions all come from package mvabund.
Import the environmental data - the one cleaned, prepared and saved in the previous notebook (classical multivariate methods). It contains long-term averages for the water column data (as long-term as available, at least) at each station, repeated for each replicate, the sediment data (2013-2014), and the seagrass data (2013-2014), again repeated to the same number of replicates. Only the variables determined to be significant by PCA are kept.
env.zostera <- read_csv(here(save.dir, "env_data_ordinations_zostera.csv"))
## convert station to factor
(env.zostera <- env.zostera %>%
mutate(station = factor(station,
levels = c("Poda", "Otmanli", "Vromos", "Gradina", "Ropotamo")))
)
Station is a factor, the rest of the variables are numeric.
Turn the zoobenthic data (minus the all-0 taxa) into a matrix - easier for the mvabund package and methods to deal with.
## there is already one subset of filtered count data (32 x 94) - use it
zoo.mvabnd.zostera <- mvabund(zoo.abnd.flt.zostera)
First, let’s see if the groups from the latent variable model (more or less equal to the clusters from the classical ordination) are valid, and which species exhibit a response.
I’m going to try something new here - 1) loose clusters from the LVM ordination, 1 = Poda-Otmanli, 2 = Vromos, 3 = Gradina-Ropotamo. 2) stations as clusters, as I did before for the seagrass data, although I don’t believe it’s valid/justified to do so… 3) another configuration of clusters from the LVM ordination: 1 = Z1-Z2, 2 = Z3, 3 = Z4, 4 = Z5.
## construct the vectors of the clusters by hand - first, situation 1 above
lvm.clusters.zostera.1 <- rep(1:3, times = c(16, 4, 12))
(lvm.clusters.zostera.1 <- factor(lvm.clusters.zostera.1))
## again, for case 2
lvm.clusters.zostera.2 <- rep(1:5, times = c(8, 8, 4, 8, 4))
(lvm.clusters.zostera.2 <- factor(lvm.clusters.zostera.2))
## again, for case 3
lvm.clusters.zostera.3 <- rep(1:4, times = c(16, 4, 8, 4))
(lvm.clusters.zostera.3 <- factor(lvm.clusters.zostera.3))
LVM clusters - case 1 Check the model assumptions. 1. Mean-variance assumption => determines the choice of family parameter. Can be checked by plotting residuals vs fits: if little pattern - the chosen mean-variance assumption is plausible.
Another way: direct plotting (variance ~ mean), for each species within each factor level.
plot(manyglm(zoo.mvabnd.zostera ~ lvm.clusters.zostera.1, family = "negative.binomial"))
meanvar.plot(zoo.mvabnd.zostera ~ lvm.clusters.zostera.1, table = TRUE)
It’s not perfect, but it’s not too terrible either.
Everything looks more or less fine; fit the model.
glms.lvm.zostera.1 <- manyglm(zoo.mvabnd.zostera ~ lvm.clusters.zostera.1,
family = "negative.binomial")
Explore the fit (residuals, diagnostic plots, etc.).
## residuals vs fitted values
plot(glms.lvm.zostera.1)
## all traditional (g)lm diagnostic plots
plot.manyglm(glms.lvm.zostera.1, which = 1:3)
# ### source mvabund GLM plotting functions modified to use a grey palette - I just can't redo these plots on my own, the function is doing too complicated things internally to scale the x and y axes
# source(here(functions.dir, "default.plot.manyglm_grey.R"))
# source(here(functions.dir, "plot.manyglm_grey.R"))
#
# par(mfrow = c(2,2))
# lapply(1:3, function(i) plot.manyglm.grey(glms.lvm.zostera, which = i, sub.caption = ""))
# par(mfrow = c(1, 1))
I really don’t like the rainbow palette, but I would like to include these plots in my thesis results.. Will have to do something about it, just not right now.
Save the model!
write_rds(glms.lvm.zostera.1,
here(save.dir, "glms_lvm_zostera_1.RDS"))
Let’s see the model summary (NB takes a LOT of time if there are many resamplings!).
(glms.lvm.zostera.1.summary <- summary(glms.lvm.zostera.1,
test = "LR", p.uni = "adjusted",
nBoot = 999, ## limit the number of permutations if you just want to check it out
show.time = "all")
)
The factor is highly significant according to the models.
This also allows us to see which species exhibit a response to the chosen factor. The LR (likelihood ratio) statistic is used as a measure of the strength of individual taxon contributions to the observed patterns. I’ll save the summary for safekeeping, but I’ll also run an anova - to get an analysis of deviance table on the model fit (also better for extracting the species contributions, or at least I know how to do it).
write_rds(glms.lvm.zostera.1.summary,
here(save.dir, "glms_lvm_zostera_1_summary.RDS"))
Run the anova on the model.
(glms.lvm.zostera.1.aov <- anova.manyglm(glms.lvm.zostera.1,
test = "LR", p.uni = "adjusted",
nBoot = 999, ## limit the number of permutations for a shorter run time
pairwise.comp = ~lvm.clusters.zostera.1, ## check the pairwise comparison between clusters
show.time = "all")
)
Resampling begins for test 1.
Resampling run 0 finished. Time elapsed: 0.00 minutes...
Resampling run 100 finished. Time elapsed: 0.12 minutes...
Resampling run 200 finished. Time elapsed: 0.25 minutes...
Resampling run 300 finished. Time elapsed: 0.38 minutes...
Resampling run 400 finished. Time elapsed: 0.50 minutes...
Resampling run 500 finished. Time elapsed: 0.63 minutes...
Resampling run 600 finished. Time elapsed: 0.75 minutes...
Resampling run 700 finished. Time elapsed: 0.88 minutes...
Resampling run 800 finished. Time elapsed: 1.01 minutes...
Resampling run 900 finished. Time elapsed: 1.13 minutes...
Time elapsed: 0 hr 1 min 14 sec
Analysis of Deviance Table
Model: manyglm(formula = zoo.mvabnd.zostera ~ lvm.clusters.zostera.1,
Model: family = "negative.binomial")
Multivariate test:
Res.Df Df.diff Dev Pr(>Dev)
(Intercept) 31
lvm.clusters.zostera.1 29 2 703.6 0.001 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Pairwise comparison results:
Observed statistic
lvm.clusters.zostera.1:1 vs lvm.clusters.zostera.1:3 370.0
lvm.clusters.zostera.1:2 vs lvm.clusters.zostera.1:3 346.1
lvm.clusters.zostera.1:1 vs lvm.clusters.zostera.1:2 319.4
Free Stepdown Adjusted P-Value
lvm.clusters.zostera.1:1 vs lvm.clusters.zostera.1:3 0.001 ***
lvm.clusters.zostera.1:2 vs lvm.clusters.zostera.1:3 0.001 ***
lvm.clusters.zostera.1:1 vs lvm.clusters.zostera.1:2 0.001 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Univariate Tests:
Abra.alba Abra.sp. Actiniaria
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 25.931 0.002 4.159 0.986 1.386 1.000
Alitta.succinea Ampelisca.diadema
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 8.646 0.434 20.059 0.005
Amphibalanus.improvisus Ampithoe.sp.
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 14.882 0.036 6.321 0.802
Anadara.kagoshimensis Apherusa.bispinosa
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 1.386 1.000 1.962 1.000
Apseudopsis.ostroumovi Bittium.reticulatum
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 30.855 0.001 31.183 0.001
Brachynotus.sexdentatus Capitella.capitata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 0.575 1.000 26.213 0.002
Capitella.minima Chamelea.gallina
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 0.739 1.000 26.23 0.002
Chironomidae.larvae Cumella.limicola
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 4.195 0.979 32.545 0.001
Cumella.pygmaea Cytharella.costulata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 4.228 0.979 3.122 0.996
Diogenes.pugilator Eteone.flava Eunice.vittata
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.1 0.896 1.000 4.159 0.987 0.871
Eurydice.dollfusi Exogone.naidina
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 1.000 4.228 0.979 5.944 0.839
Gastrosaccus.sanctus Genetyllis.tuberculata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 1.962 1.000 0.467 1.000
Glycera.sp. Glycera.tridactyla
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 7.992 0.559 4.484 0.961
Glycera.unicornis Harmothoe.imbricata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 1.386 1.000 1.841 1.000
Harmothoe.reticulata Heteromastus.filiformis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 5.13 0.925 10.222 0.247
Hirudinea Hydrobia.acuta Hydrobia.sp.
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 0.915 1.000 1.438 1.000 1.681 1.000
Iphinoe.tenella Kellia.suborbicularis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 4.16 0.981 4.584 0.949
Lagis.koreni Leiochone.leiopygos
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 9.705 0.296 6.218 0.812
Lentidium.mediterraneum Lepidochitona.cinerea
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 3.923 0.989 1.962 1.000
Loripes.orbiculatus Lucinella.divaricata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 25.259 0.002 4.257 0.979
Magelona.papillicornis Maldane.glebifex
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 1.962 1.000 1.386 1.000
Melinna.palmata Microdeutopus.gryllotalpa
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 24.203 0.003 3.088 0.996
Micromaldane.ornithochaeta Micronephthys.stammeri
Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.1 4.544 0.959 0.575
Microphthalmus.fragilis Microphthalmus.sp.
Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.1 1.000 1.437 1.000 9.396
Monocorophium.acherusicum Mytilaster.lineatus
Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.1 0.342 38.329 0.001 21.933
Mytilus.galloprovincialis Nemertea
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 0.003 1.386 1.000 3.554 0.995
Nephtys.cirrosa Nephtys.kersivalensis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 0.575 1.000 1.386 1.000
Nereis.perivisceralis Nereis.pulsatoria
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 1.962 1.000 2.987 0.998
Nototropis.guttatus Oligochaeta
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 2.61 0.999 21.776 0.003
Paradoneis.harpagonea Parthenina.interstincta
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 1.386 1.000 1.965 1.000
Parvicardium.exiguum Perinereis.cultrifera
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 2.086 1.000 3.68 0.995
Perioculodes.longimanus Phoronida
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 2.254 0.999 4.088 0.989
Phyllodoce.sp. Platyhelminthes
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 4.159 0.987 4.534 0.959
Platynereis.dumerilii Polititapes.aureus
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 7.28 0.627 3.315 0.996
Polychaeta.larvae Polydora.ciliata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 4.159 0.984 15.69 0.022
Polygordius.neapolitanus Prionospio.cirrifera
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 4.159 0.984 13.193 0.072
Protodorvillea.kefersteini Pseudocuma.longicorne
Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.1 17.536 0.009 1.491
Rissoa.membranacea Rissoa.splendida
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 1.000 0.158 1.000 11.549 0.145
Salvatoria.clavata Schistomeringos.rudolphi
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 17.003 0.015 2.634 0.999
Sphaerosyllis.hystrix Spio.filicornis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 12.713 0.085 28.109 0.001
Spisula.subtruncata Stenosoma.capito
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 5.885 0.842 5.011 0.930
Syllis.gracilis Syllis.hyalina Tellina.tenuis
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.1 5.956 0.839 2.048 1.000 1.491
Thracia.phaseolina Tricolia.pullus
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 1.000 1.962 1.000 7.956 0.559
Tritia.neritea Tritia.reticulata Turbellaria
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.1 1.406 1.000 2.932 0.999 0.135
Upogebia.pusilla
Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 1.000 4.257 0.979
Arguments:
Test statistics calculated assuming uncorrelated response (for faster computation)
P-value calculated using 999 resampling iterations via PIT-trap resampling (to account for correlation in testing.
Aaaand the differences between clusters are highly significant in this iteration..
Save the ANOVA, too.
write_rds(glms.lvm.zostera.1.aov,
here(save.dir, "glms_lvm_zostera_1_anova.RDS"))
NOW let’s get the taxa with the highest contributions to the tested pattern.
## get the top contributing species for the initial zostera GLMs
(top.sp.glms.lvm.zostera.1 <- top_n_sp_glm(glms.lvm.zostera.1.aov, tot.dev.expl = 0.75)
)
[1] "Total deviance explained: 0.76"
Monocorophium.acherusicum Cumella.limicola Bittium.reticulatum
38.328670 32.545296 31.182981
Apseudopsis.ostroumovi Spio.filicornis Chamelea.gallina
30.854772 28.108956 26.230014
Capitella.capitata Abra.alba Loripes.orbiculatus
26.212753 25.931365 25.258931
Melinna.palmata Mytilaster.lineatus Oligochaeta
24.202974 21.932810 21.776133
Ampelisca.diadema Protodorvillea.kefersteini Salvatoria.clavata
20.058951 17.536359 17.003452
Polydora.ciliata Amphibalanus.improvisus Prionospio.cirrifera
15.689516 14.882244 13.192537
Sphaerosyllis.hystrix Rissoa.splendida Heteromastus.filiformis
12.713470 11.548802 10.221754
Lagis.koreni Microphthalmus.sp. Alitta.succinea
9.704792 9.395618 8.646146
Glycera.sp. Tricolia.pullus Platynereis.dumerilii
7.992056 7.955930 7.279560
Ampithoe.sp. Leiochone.leiopygos Syllis.gracilis
6.321173 6.218039 5.956444
## unfortunately, mvabund likes to rename my species when converting the data to matrix (no spaces in names), and since I'm going to look them up in my initial untransformed count data, I have to change them back..
names(top.sp.glms.lvm.zostera.1) <- names(top.sp.glms.lvm.zostera.1) %>%
str_replace(pattern = "\\.", replacement = " ")
top.sp.glms.lvm.zostera.1
Monocorophium acherusicum Cumella limicola Bittium reticulatum
38.328670 32.545296 31.182981
Apseudopsis ostroumovi Spio filicornis Chamelea gallina
30.854772 28.108956 26.230014
Capitella capitata Abra alba Loripes orbiculatus
26.212753 25.931365 25.258931
Melinna palmata Mytilaster lineatus Oligochaeta
24.202974 21.932810 21.776133
Ampelisca diadema Protodorvillea kefersteini Salvatoria clavata
20.058951 17.536359 17.003452
Polydora ciliata Amphibalanus improvisus Prionospio cirrifera
15.689516 14.882244 13.192537
Sphaerosyllis hystrix Rissoa splendida Heteromastus filiformis
12.713470 11.548802 10.221754
Lagis koreni Microphthalmus sp. Alitta succinea
9.704792 9.395618 8.646146
Glycera sp. Tricolia pullus Platynereis dumerilii
7.992056 7.955930 7.279560
Ampithoe sp. Leiochone leiopygos Syllis gracilis
6.321173 6.218039 5.956444
Try to plot these top contributing species - for whatever that’s worth, because 50 species on a plot is still a monstrosity.
## get the species and their abundances from the original count data, and transform them to long format
(abnd.top.sp.glms.lvm.zostera.1 <- zoo.abnd.zostera %>%
select(station, names(top.sp.glms.lvm.zostera.1)) %>%
gather(key = "species", value = "count", -station) %>%
## turn species into a factor, or you'll be very very sorry later, when they're out of order on the plot. NB need to be in REVERSE order, because ggplot plots from bottom to top, and I want the top-contributing species on top.
mutate(species = factor(species, levels = rev(names(top.sp.glms.lvm.zostera.1)))) %>%
mutate(group = factor(case_when(station %in% c("Poda", "Otmanli") ~ 1,
station == "Vromos" ~ 2,
station %in% c("Gradina", "Ropotamo") ~ 3))) ## add the groups to the long df
)
(plot.top.sp.glms.lvm.zostera.1 <- plot_top_n(abnd.top.sp.glms.lvm.zostera.1,
mapping = aes(x = species, y = log_y_min(count), colour = group),
labs.legend = paste0("group", as.character(levels(abnd.top.sp.glms.lvm.zostera.1$group))),
lab.y = "Abundance (log(y/min + 1))",
palette = "Set2"
) +
theme(legend.position = "top", legend.title = element_blank())
)
Well this is a nightmarish plot, but more tolerable than the one for the sand stations - there are less species here, so at least it’s readable..
Extract the top-contributing species to each cluster (this same nightmare above, but as a table). This chunk is STILL hopelessly ugly and clumsy.
top.sp.abnd.glms.lvm.zostera.1 <- lapply(names(glms.lvm.zostera.1.summary$aliased), function(x) top_sp_glms_table(glms.lvm.zostera.1.summary, x, p = 0.05))
## fix species names (remove dot)
top.sp.abnd.glms.lvm.zostera.1 <- lapply(top.sp.abnd.glms.lvm.zostera.1, function(x) x %>% mutate(species = str_replace(species, pattern = "\\.", replacement = " ")))
## rename columns (= group names) - right now they are something like "lvm.clusters.zostera2" etc.
top.sp.abnd.glms.lvm.zostera.1 <- lapply(top.sp.abnd.glms.lvm.zostera.1, function(x) x %>% rename_at(vars(contains("lvm.clusters.zostera.1")), list(~str_replace_all(., pattern = "lvm.clusters.zostera.1", "group_"))))
top.sp.abnd.glms.lvm.zostera.1 <- lapply(top.sp.abnd.glms.lvm.zostera.1, function(x) x %>% rename_at(vars(contains("Intercept")), list(~str_replace_all(., pattern = "\\(Intercept\\)", "group_1"))))
## pull the abundances from the original count df and add to the summary glm tables
## make a long df of abundances & add clusters
zoo.abnd.zostera.long.1 <- zoo.abnd.zostera %>%
select(-c(month:replicate)) %>%
gather(key = "species", value = "count", -station) %>%
mutate(group = case_when(station %in% c("Poda", "Otmanli") ~ 1,
station == "Vromos" ~ 2,
station %in% c("Gradina", "Ropotamo") ~ 3)
)
## sum sp abundances by group; nest by group
zoo.abnd.zostera.long.1.smry <- zoo.abnd.zostera.long.1 %>%
group_by(species, group) %>%
summarise(total_count = sum(count)) %>%
group_by(group) %>%
nest()
## add the counts to the group dfs - wow that's an ugly, ugly hack. Wish I had more time to write this up properly..
top.sp.abnd.glms.lvm.zostera.1 <- map2(top.sp.abnd.glms.lvm.zostera.1, zoo.abnd.zostera.long.1.smry %>% pull(group), ~left_join(.x, zoo.abnd.zostera.long.1.smry %>% filter(group == .y) %>% unnest(), by = "species"))
## since these are sum counts over all the replicates (that's why the monstrous numbers), average them to be mean counts per group. NB different groups consist of different numbers of replicates, b.c. some groups consist of more than one station
(top.sp.abnd.glms.lvm.zostera.1 <- map2(top.sp.abnd.glms.lvm.zostera.1, c(16, 4, 12), function(x, y) x %>% mutate(mean_count = total_count/y))
)
[[1]]
[[2]]
[[3]]
NA
In this case, the model shows which species exhibit a reaction based on the chosen groups - in other words, which species are more likely to be more/less abundant in each group.
I have to say, in the case of the seagrasses and case 1 clusters, there are much fewer species that exhibit a significant response - around 10 for each group.
The LRs are lower for groups 2 and 3 - not sure if this means anything, but for group 1 they are much much higher..
For group 1 (= Z1-Z2), the species/taxa with significantly higher abundance are: Bittium reticulatum, Capitella minima, Oligochaeta, H. filiformis, Polydora ciliata, Prionospio cirrifera, R. membranacea, A. alba, A.diadema, M. acherusicum; and the only one with a significantly lower abundance - Chamelea gallina.
For group 2 (= Z3), the species with higher abundance are: M. acherusicum, S. filicornis, A.dadema. The species with lower abundance are: B. reticulatum, A. alba, Oligochaeta, S. clavata, P. ciliata, P. cirrifera, H. filiformis.
For group 3 (= Z4-Z5), the species with higher abundance are: Cumella limicola, Apseudopsis ostroumovi, Capitella capitata, Mytilaster lineatus, Loripes orbiculatus; less so, but still present - C. gallina, S. clavata. The species with lower abundance are: Abra alba, Melinna palmata (totally absent).
I’ll test each station as its own group, too (as I did before, with the classical multivariate methods) - I’m not sure how much I can trust this grouping (in particular group 3 is a bit far-fetched, if you ask me..).
LVM clusters - case 2 Check the model assumptions.
plot(manyglm(zoo.mvabnd.zostera ~ lvm.clusters.zostera.2, family = "negative.binomial"))
meanvar.plot(zoo.mvabnd.zostera ~ lvm.clusters.zostera.2, table = TRUE)
It’s not perfect, but it’s not too terrible either. I think it’s a little worse than the case 1 fit.
Everything looks more or less fine; fit the model.
glms.lvm.zostera.2 <- manyglm(zoo.mvabnd.zostera ~ lvm.clusters.zostera.2,
family = "negative.binomial")
Explore the fit (residuals, diagnostic plots, etc.).
## residuals vs fitted values
plot(glms.lvm.zostera.2)
## all traditional (g)lm diagnostic plots
plot.manyglm(glms.lvm.zostera.2, which = 1:3)
# ### source mvabund GLM plotting functions modified to use a grey palette - I just can't redo these plots on my own, the function is doing too complicated things internally to scale the x and y axes
# source(here(functions.dir, "default.plot.manyglm_grey.R"))
# source(here(functions.dir, "plot.manyglm_grey.R"))
#
# par(mfrow = c(2,2))
# lapply(2:3, function(i) plot.manyglm.grey(glms.lvm.zostera, which = i, sub.caption = ""))
# par(mfrow = c(2, 2))
Save the model!
write_rds(glms.lvm.zostera.2,
here(save.dir, "glms_lvm_zostera_2.RDS"))
Let’s see the model summary (NB takes a LOT of time if there are many resamplings!).
(glms.lvm.zostera.2.summary <- summary(glms.lvm.zostera.2,
test = "LR", p.uni = "adjusted",
nBoot = 999, ## limit the number of permutations if you just want to check it out
show.time = "all")
)
The factor is highly significant according to the models.
Again, save the summary for safekeeping, but also run an anova.
write_rds(glms.lvm.zostera.2.summary,
here(save.dir, "glms_lvm_zostera_2_summary.RDS"))
Run the anova on the model.
(glms.lvm.zostera.2.aov <- anova.manyglm(glms.lvm.zostera.2,
test = "LR", p.uni = "adjusted",
nBoot = 999, ## limit the number of permutations for a shorter run time
pairwise.comp = ~lvm.clusters.zostera.2, ## check the pairwise comparison between clusters
show.time = "all")
)
Resampling begins for test 1.
Resampling run 0 finished. Time elapsed: 0.00 minutes...
Resampling run 100 finished. Time elapsed: 0.12 minutes...
Resampling run 200 finished. Time elapsed: 0.24 minutes...
Resampling run 300 finished. Time elapsed: 0.37 minutes...
Resampling run 400 finished. Time elapsed: 0.48 minutes...
Resampling run 500 finished. Time elapsed: 0.60 minutes...
Resampling run 600 finished. Time elapsed: 0.72 minutes...
Resampling run 700 finished. Time elapsed: 0.84 minutes...
Resampling run 800 finished. Time elapsed: 0.96 minutes...
Resampling run 900 finished. Time elapsed: 1.08 minutes...
Time elapsed: 0 hr 1 min 12 sec
Analysis of Deviance Table
Model: manyglm(formula = zoo.mvabnd.zostera ~ lvm.clusters.zostera.2,
Model: family = "negative.binomial")
Multivariate test:
Res.Df Df.diff Dev Pr(>Dev)
(Intercept) 31
lvm.clusters.zostera.2 27 4 1223 0.001 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Pairwise comparison results:
Observed statistic
lvm.clusters.zostera.2:1 vs lvm.clusters.zostera.2:4 418.6
lvm.clusters.zostera.2:3 vs lvm.clusters.zostera.2:4 368.1
lvm.clusters.zostera.2:3 vs lvm.clusters.zostera.2:5 336.6
lvm.clusters.zostera.2:2 vs lvm.clusters.zostera.2:3 325.3
lvm.clusters.zostera.2:1 vs lvm.clusters.zostera.2:5 310.4
lvm.clusters.zostera.2:4 vs lvm.clusters.zostera.2:5 295.1
lvm.clusters.zostera.2:2 vs lvm.clusters.zostera.2:4 278.8
lvm.clusters.zostera.2:1 vs lvm.clusters.zostera.2:3 278.8
lvm.clusters.zostera.2:2 vs lvm.clusters.zostera.2:5 271.7
lvm.clusters.zostera.2:1 vs lvm.clusters.zostera.2:2 219.8
Free Stepdown Adjusted P-Value
lvm.clusters.zostera.2:1 vs lvm.clusters.zostera.2:4 0.002 **
lvm.clusters.zostera.2:3 vs lvm.clusters.zostera.2:4 0.008 **
lvm.clusters.zostera.2:3 vs lvm.clusters.zostera.2:5 0.009 **
lvm.clusters.zostera.2:2 vs lvm.clusters.zostera.2:3 0.009 **
lvm.clusters.zostera.2:1 vs lvm.clusters.zostera.2:5 0.009 **
lvm.clusters.zostera.2:4 vs lvm.clusters.zostera.2:5 0.009 **
lvm.clusters.zostera.2:2 vs lvm.clusters.zostera.2:4 0.009 **
lvm.clusters.zostera.2:1 vs lvm.clusters.zostera.2:3 0.009 **
lvm.clusters.zostera.2:2 vs lvm.clusters.zostera.2:5 0.009 **
lvm.clusters.zostera.2:1 vs lvm.clusters.zostera.2:2 0.009 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Univariate Tests:
Abra.alba Abra.sp. Actiniaria
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 27.461 0.002 4.159 1.000 2.773 1.000
Alitta.succinea Ampelisca.diadema
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 36.585 0.001 30.57 0.001
Amphibalanus.improvisus Ampithoe.sp.
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 23.408 0.009 12.209 0.448
Anadara.kagoshimensis Apherusa.bispinosa
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 2.773 1.000 2.773 1.000
Apseudopsis.ostroumovi Bittium.reticulatum
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 52.794 0.001 34.193 0.001
Brachynotus.sexdentatus Capitella.capitata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 2.773 1.000 33.199 0.001
Capitella.minima Chamelea.gallina
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 7.227 0.941 28.525 0.001
Chironomidae.larvae Cumella.limicola
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 9.651 0.761 41.328 0.001
Cumella.pygmaea Cytharella.costulata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 6.118 0.976 7.354 0.936
Diogenes.pugilator Eteone.flava Eunice.vittata
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.2 3.059 1.000 4.159 1.000 3.779
Eurydice.dollfusi Exogone.naidina
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 1.000 6.118 0.976 8.477 0.879
Gastrosaccus.sanctus Genetyllis.tuberculata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 2.773 1.000 7.712 0.926
Glycera.sp. Glycera.tridactyla
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 7.992 0.918 5.516 0.994
Glycera.unicornis Harmothoe.imbricata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 2.773 1.000 8.753 0.861
Harmothoe.reticulata Heteromastus.filiformis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 15.21 0.180 28.142 0.001
Hirudinea Hydrobia.acuta Hydrobia.sp.
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 1.726 1.000 2.932 1.000 5.959 0.980
Iphinoe.tenella Kellia.suborbicularis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 9.824 0.747 9.077 0.846
Lagis.koreni Leiochone.leiopygos
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 16.471 0.131 15.331 0.179
Lentidium.mediterraneum Lepidochitona.cinerea
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 5.545 0.990 2.773 1.000
Loripes.orbiculatus Lucinella.divaricata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 31.348 0.001 6.172 0.974
Magelona.papillicornis Maldane.glebifex
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 2.773 1.000 2.773 1.000
Melinna.palmata Microdeutopus.gryllotalpa
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 26.86 0.002 16.35 0.131
Micromaldane.ornithochaeta Micronephthys.stammeri
Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.2 6.655 0.962 4.159
Microphthalmus.fragilis Microphthalmus.sp.
Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.2 0.999 2.931 1.000 11.855
Monocorophium.acherusicum Mytilaster.lineatus
Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.2 0.487 39.891 0.001 47.863
Mytilus.galloprovincialis Nemertea
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 0.001 2.773 1.000 10.046 0.723
Nephtys.cirrosa Nephtys.kersivalensis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 2.773 1.000 2.773 1.000
Nereis.perivisceralis Nereis.pulsatoria
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 2.773 1.000 6.238 0.971
Nototropis.guttatus Oligochaeta
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 7.21 0.941 36.765 0.001
Paradoneis.harpagonea Parthenina.interstincta
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 2.773 1.000 4.546 0.997
Parvicardium.exiguum Perinereis.cultrifera
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 2.973 1.000 6.689 0.962
Perioculodes.longimanus Phoronida
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 4.365 0.999 5.483 0.994
Phyllodoce.sp. Platyhelminthes
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 4.159 1.000 8.508 0.879
Platynereis.dumerilii Polititapes.aureus
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 8.91 0.851 10.998 0.575
Polychaeta.larvae Polydora.ciliata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 4.159 0.999 41.905 0.001
Polygordius.neapolitanus Prionospio.cirrifera
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 4.159 0.999 46.273 0.001
Protodorvillea.kefersteini Pseudocuma.longicorne
Dev Pr(>Dev) Dev
(Intercept) <NA> <NA> <NA>
lvm.clusters.zostera.2 30.14 0.001 3.112
Rissoa.membranacea Rissoa.splendida
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA> <NA>
lvm.clusters.zostera.2 1.000 2.72 1.000 17.118 0.113
Salvatoria.clavata Schistomeringos.rudolphi
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA>
lvm.clusters.zostera.2 32.959 0.001 13.822 0.291
Sphaerosyllis.hystrix Spio.filicornis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA>
lvm.clusters.zostera.2 16.76 0.124 34.879 0.001
Spisula.subtruncata Stenosoma.capito
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA>
lvm.clusters.zostera.2 5.885 0.980 13.85 0.291
Syllis.gracilis Syllis.hyalina Tellina.tenuis
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept) <NA> <NA> <NA> <NA> <NA>
lvm.clusters.zostera.2 28.543 0.001 4.555 0.997 4.499
Thracia.phaseolina Tricolia.pullus
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA> <NA>
lvm.clusters.zostera.2 0.997 2.773 1.000 14.668 0.223
Tritia.neritea Tritia.reticulata Turbellaria
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept) <NA> <NA> <NA> <NA> <NA>
lvm.clusters.zostera.2 5.516 0.994 11.233 0.558 2.652
Upogebia.pusilla
Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA>
lvm.clusters.zostera.2 1.000 10.026 0.723
Arguments:
Test statistics calculated assuming uncorrelated response (for faster computation)
P-value calculated using 999 resampling iterations via PIT-trap resampling (to account for correlation in testing.
Again, these groups are sufficiently different from one another.. No clue here.
Save the ANOVA, too.
write_rds(glms.lvm.zostera.2.aov,
here(save.dir, "glms_lvm_zostera_2_anova.RDS"))
NOW let’s get the taxa with the highest contributions to the tested pattern.
## get the top contributing species for the initial zostera GLMs
(top.sp.glms.lvm.zostera.2 <- top_n_sp_glm(glms.lvm.zostera.2.aov, tot.dev.expl = 0.75)
)
[1] "Total deviance explained: 0.799"
Apseudopsis.ostroumovi Mytilaster.lineatus Prionospio.cirrifera
52.793744 47.863335 46.272692
Polydora.ciliata Cumella.limicola Monocorophium.acherusicum
41.905219 41.327723 39.891356
Oligochaeta Alitta.succinea Spio.filicornis
36.764590 36.584687 34.878811
Bittium.reticulatum Capitella.capitata Salvatoria.clavata
34.192799 33.198554 32.958829
Loripes.orbiculatus Ampelisca.diadema Protodorvillea.kefersteini
31.347733 30.569510 30.139960
Syllis.gracilis Chamelea.gallina Heteromastus.filiformis
28.542915 28.524575 28.142266
Abra.alba Melinna.palmata Amphibalanus.improvisus
27.461185 26.859546 23.408099
Rissoa.splendida Sphaerosyllis.hystrix Lagis.koreni
17.117581 16.760313 16.471224
Microdeutopus.gryllotalpa Leiochone.leiopygos Harmothoe.reticulata
16.349670 15.331303 15.209809
Tricolia.pullus Stenosoma.capito Schistomeringos.rudolphi
14.667741 13.849706 13.822279
Ampithoe.sp. Microphthalmus.sp. Tritia.reticulata
12.208710 11.854688 11.232600
Polititapes.aureus Nemertea Upogebia.pusilla
10.998455 10.045834 10.025894
Iphinoe.tenella Chironomidae.larvae Kellia.suborbicularis
9.824385 9.650905 9.077281
Platynereis.dumerilii
8.910425
## unfortunately, mvabund likes to rename my species when converting the data to matrix (no spaces in names), and since I'm going to look them up in my initial untransformed count data, I have to change them back..
names(top.sp.glms.lvm.zostera.2) <- names(top.sp.glms.lvm.zostera.2) %>%
str_replace(pattern = "\\.", replacement = " ")
top.sp.glms.lvm.zostera.2
Apseudopsis ostroumovi Mytilaster lineatus Prionospio cirrifera
52.793744 47.863335 46.272692
Polydora ciliata Cumella limicola Monocorophium acherusicum
41.905219 41.327723 39.891356
Oligochaeta Alitta succinea Spio filicornis
36.764590 36.584687 34.878811
Bittium reticulatum Capitella capitata Salvatoria clavata
34.192799 33.198554 32.958829
Loripes orbiculatus Ampelisca diadema Protodorvillea kefersteini
31.347733 30.569510 30.139960
Syllis gracilis Chamelea gallina Heteromastus filiformis
28.542915 28.524575 28.142266
Abra alba Melinna palmata Amphibalanus improvisus
27.461185 26.859546 23.408099
Rissoa splendida Sphaerosyllis hystrix Lagis koreni
17.117581 16.760313 16.471224
Microdeutopus gryllotalpa Leiochone leiopygos Harmothoe reticulata
16.349670 15.331303 15.209809
Tricolia pullus Stenosoma capito Schistomeringos rudolphi
14.667741 13.849706 13.822279
Ampithoe sp. Microphthalmus sp. Tritia reticulata
12.208710 11.854688 11.232600
Polititapes aureus Nemertea Upogebia pusilla
10.998455 10.045834 10.025894
Iphinoe tenella Chironomidae larvae Kellia suborbicularis
9.824385 9.650905 9.077281
Platynereis dumerilii
8.910425
Try to plot these top contributing species - for whatever that’s worth, because 50 species on a plot is still a monstrosity.
## get the species and their abundances from the original count data, and transform them to long format
(abnd.top.sp.glms.lvm.zostera.2 <- zoo.abnd.zostera %>%
select(station, names(top.sp.glms.lvm.zostera.2)) %>%
gather(key = "species", value = "count", -station) %>%
## turn species into a factor, or you'll be very very sorry later, when they're out of order on the plot. NB need to be in REVERSE order, because ggplot plots from bottom to top, and I want the top-contributing species on top.
mutate(species = factor(species, levels = rev(names(top.sp.glms.lvm.zostera.2)))) %>%
mutate(group = factor(case_when(station == "Poda" ~ 1,
station == "Otmanli" ~ 2,
station == "Vromos" ~ 3,
station == "Gradina" ~ 4,
station == "Ropotamo" ~ 5)))
)
(plot.top.sp.glms.lvm.zostera.2 <- plot_top_n(abnd.top.sp.glms.lvm.zostera.2,
mapping = aes(x = species, y = log_y_min(count), colour = group),
labs.legend = paste0("group", as.character(levels(abnd.top.sp.glms.lvm.zostera.2$group))),
lab.y = "Abundance (log(y/min + 1))",
palette = "Set2"
) +
theme(legend.position = "top", legend.title = element_blank())
)
Extract the top-contributing species to each cluster (this same nightmare above, but as a table). This chunk is STILL hopelessly ugly and clumsy.
top.sp.abnd.glms.lvm.zostera.2 <- lapply(names(glms.lvm.zostera.2.summary$aliased), function(x) top_sp_glms_table(glms.lvm.zostera.2.summary, x, p = 0.05))
## fix species names (remove dot)
top.sp.abnd.glms.lvm.zostera.2 <- lapply(top.sp.abnd.glms.lvm.zostera.2, function(x) x %>% mutate(species = str_replace(species, pattern = "\\.", replacement = " ")))
## rename columns (= group names) - right now they are something like "lvm.clusters.zostera2" etc.
top.sp.abnd.glms.lvm.zostera.2 <- lapply(top.sp.abnd.glms.lvm.zostera.2, function(x) x %>% rename_at(vars(contains("lvm.clusters.zostera.2")), list(~str_replace_all(., pattern = "lvm.clusters.zostera.2", "group_"))))
top.sp.abnd.glms.lvm.zostera.2 <- lapply(top.sp.abnd.glms.lvm.zostera.2, function(x) x %>% rename_at(vars(contains("Intercept")), list(~str_replace_all(., pattern = "\\(Intercept\\)", "group_1"))))
## pull the abundances from the original count df and add to the summary glm tables
## make a long df of abundances & add clusters
zoo.abnd.zostera.long.2 <- zoo.abnd.zostera %>%
select(-c(month:replicate)) %>%
gather(key = "species", value = "count", -station) %>%
mutate(group = case_when(station == "Poda" ~ 1,
station == "Otmanli" ~ 2,
station == "Vromos" ~ 3,
station == "Gradina" ~ 4,
station == "Ropotamo" ~ 5)
)
## sum sp abundances by group; nest by group
zoo.abnd.zostera.long.2.smry <- zoo.abnd.zostera.long.2 %>%
group_by(species, group) %>%
summarise(total_count = sum(count)) %>%
group_by(group) %>%
nest()
## add the counts to the group dfs - wow that's an ugly, ugly hack. Wish I had more time to write this up properly..
top.sp.abnd.glms.lvm.zostera.2 <- map2(top.sp.abnd.glms.lvm.zostera.2, zoo.abnd.zostera.long.2.smry %>% pull(group), ~left_join(.x, zoo.abnd.zostera.long.2.smry %>% filter(group == .y) %>% unnest(), by = "species"))
## since these are sum counts over all the replicates (that's why the monstrous numbers), average them to be mean counts per group. NB different groups consist of different numbers of replicates, b.c. some groups consist of more than one station
(top.sp.abnd.glms.lvm.zostera.2 <- map2(top.sp.abnd.glms.lvm.zostera.2, c(8, 8, 4, 8, 4), function(x, y) x %>% mutate(mean_count = total_count/y))
)
[[1]]
[[2]]
[[3]]
[[4]]
[[5]]
NA
In the case of the seagrasses and case 2 clusters (= stations), the picture is still more unclear.. I suppose this is in no small part because of the differences 2013-14 - very marked for Poda and Otmanli. I suspect the stations changed in these two years (we were looking for Z. noltii in 2014 in particular) - but still, there is much variability. In the future, it’s probably going to be worth it to have more stations in a meadow, if we really want to have an idea of the communities there, and their variability.
The LRs seem to be a bit lower for groups 2, 4, maybe 5 too - still not sure if you can use that as a significance measure.
For now, in group 1 (= Z1), it’s hard to pick some characteristic species - because of the variability between 2013-2014, no doubt. The species/taxa with significantly higher abundance are: Bittium reticulatum, Capitella minima, Polydora ciliata, Prionosprio cirrifera (+ others, medium abundance); and the ones with a significantly lower abundance - or even absent - C. gallina, A. ostroumovi, S. clavata, C. limicola, C. costulata, S. hystrix, S. gracilis, T. pullus.
For group 2 (= Z2), the species with higher abundance - which is not really all that high; this group is also loose, hard to distinguish from group 1 - are: S. gracilis, M. lineatus, P. ciliata. The only species with lower abundance - in fact 0 - is Alitta succinea.
For group 3 (= Z3), the species with higher abundance are: M. acherusicum, S. filicornis, A. diadema. The species with lower abundance (or 0) are: B. reticulatum, P. ciliata, P. cirrifera, A. alba, A. succiena, S. clavata, Oligochaeta, A. improvisus.
For group 4 (= Z4), the species with higher abundance are: M. lineatus (very much so); C. limicola, P. kefersteini, C. gallina, C. capitata. The species with lower abundance (or 0) are: P. ciliata, P. cirrifera, A. succinea, A. improvisus, A. alba.
For group 5 (= Z5), the species with higher abundance are: A. ostroumovi, C. capitata, Oligochaeta. The species with lower abundance (or 0) are: R. splendida, T. pullus.
LVM clusters - case 3 Last try: group 1 = Z1-Z2, group 2 = Z3, group 3 = Z4, group 4 = Z5.
Check the model assumptions.
plot(manyglm(zoo.mvabnd.zostera ~ lvm.clusters.zostera.3, family = "negative.binomial"))
meanvar.plot(zoo.mvabnd.zostera ~ lvm.clusters.zostera.3, table = TRUE)
More or less the same as case 2 before it.
Everything looks more or less fine; fit the model.
glms.lvm.zostera.3 <- manyglm(zoo.mvabnd.zostera ~ lvm.clusters.zostera.3,
family = "negative.binomial")
Explore the fit (residuals, diagnostic plots, etc.).
## residuals vs fitted values
plot(glms.lvm.zostera.3)
## all traditional (g)lm diagnostic plots
plot.manyglm(glms.lvm.zostera.3, which = 1:3)
# ### source mvabund GLM plotting functions modified to use a grey palette - I just can't redo these plots on my own, the function is doing too complicated things internally to scale the x and y axes
# source(here(functions.dir, "default.plot.manyglm_grey.R"))
# source(here(functions.dir, "plot.manyglm_grey.R"))
#
# par(mfrow = c(3,3))
# lapply(3:3, function(i) plot.manyglm.grey(glms.lvm.zostera, which = i, sub.caption = ""))
# par(mfrow = c(3, 3))
Save the model!
write_rds(glms.lvm.zostera.3,
here(save.dir, "glms_lvm_zostera_3.RDS"))
Let’s see the model summary (NB takes a LOT of time if there are many resamplings!).
(glms.lvm.zostera.3.summary <- summary(glms.lvm.zostera.3,
test = "LR", p.uni = "adjusted",
nBoot = 999, ## limit the number of permutations if you just want to check it out
show.time = "all")
)
The factor is highly significant according to the models.
Again, save the summary for safekeeping, but also run an anova.
write_rds(glms.lvm.zostera.3.summary,
here(save.dir, "glms_lvm_zostera_3_summary.RDS"))
Run the anova on the model.
(glms.lvm.zostera.3.aov <- anova.manyglm(glms.lvm.zostera.3,
test = "LR", p.uni = "adjusted",
nBoot = 999, ## limit the number of permutations for a shorter run time
pairwise.comp = ~lvm.clusters.zostera.3, ## check the pairwise comparisons between clusters
show.time = "all")
)
Resampling begins for test 1.
Resampling run 0 finished. Time elapsed: 0.00 minutes...
Resampling run 100 finished. Time elapsed: 0.12 minutes...
Resampling run 200 finished. Time elapsed: 0.24 minutes...
Resampling run 300 finished. Time elapsed: 0.37 minutes...
Resampling run 400 finished. Time elapsed: 0.50 minutes...
Resampling run 500 finished. Time elapsed: 0.63 minutes...
Resampling run 600 finished. Time elapsed: 0.75 minutes...
Resampling run 700 finished. Time elapsed: 0.88 minutes...
Resampling run 800 finished. Time elapsed: 1.00 minutes...
Resampling run 900 finished. Time elapsed: 1.13 minutes...
Time elapsed: 0 hr 1 min 15 sec
Analysis of Deviance Table
Model: manyglm(formula = zoo.mvabnd.zostera ~ lvm.clusters.zostera.3,
Model: family = "negative.binomial")
Multivariate test:
Res.Df Df.diff Dev Pr(>Dev)
(Intercept) 31
lvm.clusters.zostera.3 28 3 965.7 0.001 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Pairwise comparison results:
Observed statistic
lvm.clusters.zostera.3:2 vs lvm.clusters.zostera.3:3 368.1
lvm.clusters.zostera.3:1 vs lvm.clusters.zostera.3:3 355.8
lvm.clusters.zostera.3:2 vs lvm.clusters.zostera.3:4 336.6
lvm.clusters.zostera.3:1 vs lvm.clusters.zostera.3:2 319.4
lvm.clusters.zostera.3:3 vs lvm.clusters.zostera.3:4 295.1
lvm.clusters.zostera.3:1 vs lvm.clusters.zostera.3:4 275.0
Free Stepdown Adjusted P-Value
lvm.clusters.zostera.3:2 vs lvm.clusters.zostera.3:3 0.007 **
lvm.clusters.zostera.3:1 vs lvm.clusters.zostera.3:3 0.007 **
lvm.clusters.zostera.3:2 vs lvm.clusters.zostera.3:4 0.010 **
lvm.clusters.zostera.3:1 vs lvm.clusters.zostera.3:2 0.010 **
lvm.clusters.zostera.3:3 vs lvm.clusters.zostera.3:4 0.010 **
lvm.clusters.zostera.3:1 vs lvm.clusters.zostera.3:4 0.010 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Univariate Tests:
Abra.alba Abra.sp. Actiniaria
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 27.438 0.002 4.159 1.000 1.386 1.000
Alitta.succinea Ampelisca.diadema
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 13.528 0.170 20.153 0.014
Amphibalanus.improvisus Ampithoe.sp.
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 15.231 0.098 8.727 0.657
Anadara.kagoshimensis Apherusa.bispinosa
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 1.386 1.000 2.773 1.000
Apseudopsis.ostroumovi Bittium.reticulatum
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 51.414 0.001 34.189 0.001
Brachynotus.sexdentatus Capitella.capitata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 1.386 1.000 30.336 0.001
Capitella.minima Chamelea.gallina
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 6.793 0.868 27.138 0.002
Chironomidae.larvae Cumella.limicola
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 9.651 0.539 33.634 0.001
Cumella.pygmaea Cytharella.costulata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 6.118 0.940 3.195 1.000
Diogenes.pugilator Eteone.flava Eunice.vittata
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.3 2.38 1.000 4.159 1.000 3.572
Eurydice.dollfusi Exogone.naidina
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 1.000 6.118 0.940 7.227 0.845
Gastrosaccus.sanctus Genetyllis.tuberculata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 2.773 1.000 1.966 1.000
Glycera.sp. Glycera.tridactyla
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 7.992 0.753 4.484 0.997
Glycera.unicornis Harmothoe.imbricata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 1.386 1.000 3.192 1.000
Harmothoe.reticulata Heteromastus.filiformis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 9.416 0.574 17.463 0.039
Hirudinea Hydrobia.acuta Hydrobia.sp.
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 1.726 1.000 1.438 1.000 2.25 1.000
Iphinoe.tenella Kellia.suborbicularis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 5.47 0.985 8.973 0.633
Lagis.koreni Leiochone.leiopygos
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 10.472 0.426 14.995 0.102
Lentidium.mediterraneum Lepidochitona.cinerea
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 5.545 0.978 2.773 1.000
Loripes.orbiculatus Lucinella.divaricata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 28.176 0.002 6.172 0.935
Magelona.papillicornis Maldane.glebifex
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 2.773 1.000 1.386 1.000
Melinna.palmata Microdeutopus.gryllotalpa
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 24.203 0.006 13.747 0.156
Micromaldane.ornithochaeta Micronephthys.stammeri
Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.3 6.655 0.888 2.773
Microphthalmus.fragilis Microphthalmus.sp.
Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.3 1.000 1.437 1.000 11.855
Monocorophium.acherusicum Mytilaster.lineatus
Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.3 0.280 39.678 0.001 34.91
Mytilus.galloprovincialis Nemertea
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 0.001 1.386 1.000 4.566 0.996
Nephtys.cirrosa Nephtys.kersivalensis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 1.386 1.000 1.386 1.000
Nereis.perivisceralis Nereis.pulsatoria
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 2.773 1.000 2.987 1.000
Nototropis.guttatus Oligochaeta
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 7.095 0.848 27.378 0.002
Paradoneis.harpagonea Parthenina.interstincta
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 1.386 1.000 2.884 1.000
Parvicardium.exiguum Perinereis.cultrifera
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 2.911 1.000 4.654 0.994
Perioculodes.longimanus Phoronida
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 4.365 0.997 5.372 0.987
Phyllodoce.sp. Platyhelminthes
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 4.159 1.000 5.906 0.946
Platynereis.dumerilii Polititapes.aureus
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 8.537 0.673 4.065 1.000
Polychaeta.larvae Polydora.ciliata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 4.159 1.000 29.489 0.002
Polygordius.neapolitanus Prionospio.cirrifera
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 4.159 1.000 41.351 0.001
Protodorvillea.kefersteini Pseudocuma.longicorne
Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.3 25.261 0.006 1.726
Rissoa.membranacea Rissoa.splendida
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 1.000 1.701 1.000 17.118 0.046
Salvatoria.clavata Schistomeringos.rudolphi
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 32.959 0.001 5.024 0.990
Sphaerosyllis.hystrix Spio.filicornis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 12.736 0.208 34.632 0.001
Spisula.subtruncata Stenosoma.capito
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 5.885 0.946 7.189 0.848
Syllis.gracilis Syllis.hyalina Tellina.tenuis
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.3 10.712 0.392 4.555 0.996 3.112
Thracia.phaseolina Tricolia.pullus
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 1.000 2.773 1.000 13.288 0.177
Tritia.neritea Tritia.reticulata Turbellaria
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept) <NA> <NA> <NA> <NA> <NA>
lvm.clusters.zostera.3 2.037 1.000 2.932 1.000 1.527
Upogebia.pusilla
Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA>
lvm.clusters.zostera.3 1.000 10.026 0.466
Arguments:
Test statistics calculated assuming uncorrelated response (for faster computation)
P-value calculated using 999 resampling iterations via PIT-trap resampling (to account for correlation in testing.
According to the pairwise comparison, the case 3 clusters are significantly different from one another.. apparently sufficiently so.
Save the ANOVA, too.
write_rds(glms.lvm.zostera.3.aov,
here(save.dir, "glms_lvm_zostera_3_anova.RDS"))
NOW let’s get the taxa with the highest contributions to the tested pattern.
## get the top contributing species for the initial zostera GLMs
(top.sp.glms.lvm.zostera.3 <- top_n_sp_glm(glms.lvm.zostera.3.aov, tot.dev.expl = 0.75)
)
[1] "Total deviance explained: 0.82"
Apseudopsis.ostroumovi Prionospio.cirrifera Monocorophium.acherusicum
51.413970 41.350974 39.678240
Mytilaster.lineatus Spio.filicornis Bittium.reticulatum
34.909870 34.632433 34.188511
Cumella.limicola Salvatoria.clavata Capitella.capitata
33.633823 32.958829 30.336457
Polydora.ciliata Loripes.orbiculatus Abra.alba
29.489490 28.175552 27.437828
Oligochaeta Chamelea.gallina Protodorvillea.kefersteini
27.378351 27.138297 25.260955
Melinna.palmata Ampelisca.diadema Heteromastus.filiformis
24.202974 20.153373 17.463258
Rissoa.splendida Amphibalanus.improvisus Leiochone.leiopygos
17.117581 15.231072 14.994968
Microdeutopus.gryllotalpa Alitta.succinea Tricolia.pullus
13.746930 13.527612 13.287785
Sphaerosyllis.hystrix Microphthalmus.sp. Syllis.gracilis
12.736427 11.854688 10.711746
Lagis.koreni Upogebia.pusilla Chironomidae.larvae
10.471907 10.025894 9.650905
Harmothoe.reticulata Kellia.suborbicularis Ampithoe.sp.
9.416160 8.973262 8.726571
Platynereis.dumerilii Glycera.sp. Exogone.naidina
8.536643 7.992056 7.227403
Stenosoma.capito Nototropis.guttatus Capitella.minima
7.189015 7.094867 6.793482
Micromaldane.ornithochaeta
6.655076
## unfortunately, mvabund likes to rename my species when converting the data to matrix (no spaces in names), and since I'm going to look them up in my initial untransformed count data, I have to change them back..
names(top.sp.glms.lvm.zostera.3) <- names(top.sp.glms.lvm.zostera.3) %>%
str_replace(pattern = "\\.", replacement = " ")
top.sp.glms.lvm.zostera.3
Apseudopsis ostroumovi Prionospio cirrifera Monocorophium acherusicum
51.413970 41.350974 39.678240
Mytilaster lineatus Spio filicornis Bittium reticulatum
34.909870 34.632433 34.188511
Cumella limicola Salvatoria clavata Capitella capitata
33.633823 32.958829 30.336457
Polydora ciliata Loripes orbiculatus Abra alba
29.489490 28.175552 27.437828
Oligochaeta Chamelea gallina Protodorvillea kefersteini
27.378351 27.138297 25.260955
Melinna palmata Ampelisca diadema Heteromastus filiformis
24.202974 20.153373 17.463258
Rissoa splendida Amphibalanus improvisus Leiochone leiopygos
17.117581 15.231072 14.994968
Microdeutopus gryllotalpa Alitta succinea Tricolia pullus
13.746930 13.527612 13.287785
Sphaerosyllis hystrix Microphthalmus sp. Syllis gracilis
12.736427 11.854688 10.711746
Lagis koreni Upogebia pusilla Chironomidae larvae
10.471907 10.025894 9.650905
Harmothoe reticulata Kellia suborbicularis Ampithoe sp.
9.416160 8.973262 8.726571
Platynereis dumerilii Glycera sp. Exogone naidina
8.536643 7.992056 7.227403
Stenosoma capito Nototropis guttatus Capitella minima
7.189015 7.094867 6.793482
Micromaldane ornithochaeta
6.655076
Try to plot these top contributing species - for whatever that’s worth, because 50 species on a plot is still a monstrosity.
## get the species and their abundances from the original count data, and transform them to long format
(abnd.top.sp.glms.lvm.zostera.3 <- zoo.abnd.zostera %>%
select(station, names(top.sp.glms.lvm.zostera.3)) %>%
gather(key = "species", value = "count", -station) %>%
## turn species into a factor, or you'll be very very sorry later, when they're out of order on the plot. NB need to be in REVERSE order, because ggplot plots from bottom to top, and I want the top-contributing species on top.
mutate(species = factor(species, levels = rev(names(top.sp.glms.lvm.zostera.3)))) %>%
mutate(group = factor(case_when(station %in% c("Poda", "Otmanli") ~ 1,
station == "Vromos" ~ 2,
station == "Gradina" ~ 3,
station == "Ropotamo" ~ 4)) ## add the groups
)
)
(plot.top.sp.glms.lvm.zostera.3 <- plot_top_n(abnd.top.sp.glms.lvm.zostera.3,
mapping = aes(x = species, y = log_y_min(count), colour = group),
labs.legend = paste0("group", as.character(levels(abnd.top.sp.glms.lvm.zostera.3$group))),
lab.y = "Abundance (log(y/min + 1))",
palette = "Set2"
) +
theme(legend.position = "top", legend.title = element_blank())
)
Extract the top-contributing species to each cluster (this same nightmare above, but as a table). This chunk is STILL hopelessly ugly and clumsy.
top.sp.abnd.glms.lvm.zostera.3 <- lapply(names(glms.lvm.zostera.3.summary$aliased), function(x) top_sp_glms_table(glms.lvm.zostera.3.summary, x, p = 0.05))
## fix species names (remove dot)
top.sp.abnd.glms.lvm.zostera.3 <- lapply(top.sp.abnd.glms.lvm.zostera.3, function(x) x %>% mutate(species = str_replace(species, pattern = "\\.", replacement = " ")))
## rename columns (= group names) - right now they are something like "lvm.clusters.zostera3" etc.
top.sp.abnd.glms.lvm.zostera.3 <- lapply(top.sp.abnd.glms.lvm.zostera.3, function(x) x %>% rename_at(vars(contains("lvm.clusters.zostera.3")), list(~str_replace_all(., pattern = "lvm.clusters.zostera.3", "group_"))))
top.sp.abnd.glms.lvm.zostera.3 <- lapply(top.sp.abnd.glms.lvm.zostera.3, function(x) x %>% rename_at(vars(contains("Intercept")), list(~str_replace_all(., pattern = "\\(Intercept\\)", "group_1"))))
## pull the abundances from the original count df and add to the summary glm tables
## make a long df of abundances & add clusters
zoo.abnd.zostera.long.3 <- zoo.abnd.zostera %>%
select(-c(month:replicate)) %>%
gather(key = "species", value = "count", -station) %>%
mutate(group = case_when(station %in% c("Poda", "Otmanli") ~ 1,
station == "Vromos" ~ 2,
station == "Gradina" ~ 3,
station == "Ropotamo" ~ 4)
)
## sum sp abundances by group; nest by group
zoo.abnd.zostera.long.3.smry <- zoo.abnd.zostera.long.3 %>%
group_by(species, group) %>%
summarise(total_count = sum(count)) %>%
group_by(group) %>%
nest()
## add the counts to the group dfs - wow that's an ugly, ugly hack. Wish I had more time to write this up properly..
top.sp.abnd.glms.lvm.zostera.3 <- map2(top.sp.abnd.glms.lvm.zostera.3, zoo.abnd.zostera.long.3.smry %>% pull(group), ~left_join(.x, zoo.abnd.zostera.long.3.smry %>% filter(group == .y) %>% unnest(), by = "species"))
## since these are sum counts over all the replicates (that's why the monstrous numbers), average them to be mean counts per group. NB different groups consist of different numbers of replicates, b.c. some groups consist of more than one station
(top.sp.abnd.glms.lvm.zostera.3 <- map2(top.sp.abnd.glms.lvm.zostera.3, c(16, 4, 8, 4), function(x, y) x %>% mutate(mean_count = total_count/y))
)
[[1]]
[[2]]
[[3]]
[[4]]
NA
In the case of the seagrasses and case 3 clusters, the picture is still more confusing..
The LRs seem to be a bit lower for groups 2 and 3, maybe 4 too - still not sure if you can use that as a significance measure.
For now, in group 1 (= Z1-Z2), the species/taxa with significantly higher abundance are: Bittium reticulatum, Capitella minima, Oligochaeta, Heteromastus filiformis, Polydora ciliata, Prionosprio cirrifera, Rissoa membranacea, Abra alba, Ampelisca diadema (+ others, medium abundance); and the ones with a significantly lower abundance - or even absent - S. clavata, A. ostroumovi, C. gallina, T. pullus.
There are more species singled out for this cluster, probably because of the variability between the two years of sampling.
For group 2 (= Z3), the species with higher abundance are: M. acherusicum, S. filicornis, A. diadema. The species with lower abundance - in fact 0 - are B. reticulatum, P. cirrifera, S. clavata, A. alba, P. ciliata, oligochaetes, H. filiformis, R. splendida.
For group 3 (= Z4), the species with higher abundance are: M. lineatus, less so - C. limicola, L. orbiculatus, P. kefersteini, C. capitata, etc. The species with lower abundance (or 0) are: P. cirrifera, P. ciliata, A. alba, etc.
For group 4 (= Z5), the species with higher abundance are: A. ostoumovi, C. capitata, oligochaetes (very abundant, but with a small LR - nice!). The species with lower abundance (or 0) are: R. splendida, S. clavata, C. limicola.
All in all, I think that group 1 (Poda + Otmanli) holds, and so does group 2 (Vromos). The question is whether to separate Gradina and Ropotamo into 2 groups, or if they make more sense together. Ropotamo is characterized by a very high number of oligochaetes, while Gradina’s most distinguishing characteristic is the high number of M. lineatus - mostly very small ones, attached to the rhizomes, close to the sediment surface I presume. Both stations have C. limicola in medium abundance, a species that is not present anywhere else (is it?).
Try to compare the three models..
(glms.lvm.zostera.comp <- anova(glms.lvm.zostera.1,
glms.lvm.zostera.2,
glms.lvm.zostera.3,
p.uni = "adjusted")
)
Time elapsed: 0 hr 2 min 22 sec
Analysis of Deviance Table
glms.lvm.zostera.1: zoo.mvabnd.zostera ~ lvm.clusters.zostera.1
glms.lvm.zostera.3: zoo.mvabnd.zostera ~ lvm.clusters.zostera.3
glms.lvm.zostera.2: zoo.mvabnd.zostera ~ lvm.clusters.zostera.2
Multivariate test:
Res.Df Df.diff Dev Pr(>Dev)
glms.lvm.zostera.1 29
glms.lvm.zostera.3 28 1 262.1 0.001 ***
glms.lvm.zostera.2 27 1 256.8 0.001 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Univariate Tests:
Abra.alba Abra.sp. Actiniaria
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 1.506 0.998 0 1.000 0 1.000
glms.lvm.zostera.2 0.023 1.000 0 1.000 1.386 1.000
Alitta.succinea Ampelisca.diadema
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 4.881 0.681 0.094 1.000
glms.lvm.zostera.2 23.057 0.001 10.416 0.104
Amphibalanus.improvisus Ampithoe.sp.
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 0.349 1.000 2.405 0.982
glms.lvm.zostera.2 8.177 0.251 3.482 0.917
Anadara.kagoshimensis Apherusa.bispinosa
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 0 1.000 0.811 1.000
glms.lvm.zostera.2 1.386 1.000 0 1.000
Apseudopsis.ostroumovi Bittium.reticulatum
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 20.559 0.001 3.006 0.941
glms.lvm.zostera.2 1.38 1.000 0.004 1.000
Brachynotus.sexdentatus Capitella.capitata
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 0.811 1.000 4.124 0.814
glms.lvm.zostera.2 1.386 1.000 2.862 0.958
Capitella.minima Chamelea.gallina
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 6.055 0.518 0.908 1.000
glms.lvm.zostera.2 0.433 1.000 1.386 0.996
Chironomidae.larvae Cumella.limicola
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 5.456 0.615 1.089 1.000
glms.lvm.zostera.2 0 1.000 7.694 0.306
Cumella.pygmaea Cytharella.costulata
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 1.891 0.996 0.073 1.000
glms.lvm.zostera.2 0 1.000 4.159 0.834
Diogenes.pugilator Eteone.flava Eunice.vittata
Dev Pr(>Dev) Dev Pr(>Dev) Dev
glms.lvm.zostera.1
glms.lvm.zostera.3 1.483 0.999 0 1.000 2.7
glms.lvm.zostera.2 0.68 1.000 0 1.000 0.208
Eurydice.dollfusi Exogone.naidina
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 0.959 1.891 0.996 1.284 1.000
glms.lvm.zostera.2 1.000 0 1.000 1.25 1.000
Gastrosaccus.sanctus Genetyllis.tuberculata
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 0.811 1.000 1.499 0.998
glms.lvm.zostera.2 0 1.000 5.746 0.574
Glycera.sp. Glycera.tridactyla Glycera.unicornis
Dev Pr(>Dev) Dev Pr(>Dev) Dev
glms.lvm.zostera.1
glms.lvm.zostera.3 0 1.000 0 1.000 0
glms.lvm.zostera.2 0 1.000 1.033 1.000 1.386
Harmothoe.imbricata Harmothoe.reticulata
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 1.000 1.352 1.000 4.286 0.796
glms.lvm.zostera.2 1.000 5.561 0.595 5.794 0.566
Heteromastus.filiformis Hirudinea Hydrobia.acuta
Dev Pr(>Dev) Dev Pr(>Dev) Dev
glms.lvm.zostera.1
glms.lvm.zostera.3 7.242 0.334 0.811 1.000 0
glms.lvm.zostera.2 10.679 0.088 0 1.000 1.494
Hydrobia.sp. Iphinoe.tenella
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 1.000 0.569 1.000 1.31 1.000
glms.lvm.zostera.2 0.996 3.709 0.906 4.354 0.811
Kellia.suborbicularis Lagis.koreni
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 4.389 0.796 0.767 1.000
glms.lvm.zostera.2 0.104 1.000 5.999 0.554
Leiochone.leiopygos Lentidium.mediterraneum
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 8.777 0.194 1.622 0.997
glms.lvm.zostera.2 0.336 1.000 0 1.000
Lepidochitona.cinerea Loripes.orbiculatus
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 0.811 1.000 2.917 0.948
glms.lvm.zostera.2 0 1.000 3.172 0.936
Lucinella.divaricata Magelona.papillicornis
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 1.915 0.995 0.811 1.000
glms.lvm.zostera.2 0 1.000 0 1.000
Maldane.glebifex Melinna.palmata
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 0 1.000 0 1.000
glms.lvm.zostera.2 1.386 0.999 2.657 0.983
Microdeutopus.gryllotalpa Micromaldane.ornithochaeta
Dev Pr(>Dev) Dev
glms.lvm.zostera.1
glms.lvm.zostera.3 10.659 0.065 2.111
glms.lvm.zostera.2 2.603 0.984 0
Micronephthys.stammeri Microphthalmus.fragilis
Pr(>Dev) Dev Pr(>Dev) Dev
glms.lvm.zostera.1
glms.lvm.zostera.3 0.995 2.197 0.994 0
glms.lvm.zostera.2 1.000 1.386 1.000 1.494
Microphthalmus.sp. Monocorophium.acherusicum
Pr(>Dev) Dev Pr(>Dev) Dev
glms.lvm.zostera.1
glms.lvm.zostera.3 1.000 2.459 0.979 1.35
glms.lvm.zostera.2 0.996 0 1.000 0.213
Mytilaster.lineatus Mytilus.galloprovincialis
Pr(>Dev) Dev Pr(>Dev) Dev
glms.lvm.zostera.1
glms.lvm.zostera.3 1.000 12.977 0.015 0
glms.lvm.zostera.2 1.000 12.953 0.038 1.386
Nemertea Nephtys.cirrosa
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 1.000 1.013 1.000 0.811 1.000
glms.lvm.zostera.2 1.000 5.479 0.634 1.386 1.000
Nephtys.kersivalensis Nereis.perivisceralis
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 0 1.000 0.811 1.000
glms.lvm.zostera.2 1.386 1.000 0 1.000
Nereis.pulsatoria Nototropis.guttatus Oligochaeta
Dev Pr(>Dev) Dev Pr(>Dev) Dev
glms.lvm.zostera.1
glms.lvm.zostera.3 0 1.000 4.485 0.751 5.602
glms.lvm.zostera.2 3.251 0.932 0.116 1.000 9.386
Paradoneis.harpagonea Parthenina.interstincta
Pr(>Dev) Dev Pr(>Dev) Dev
glms.lvm.zostera.1
glms.lvm.zostera.3 0.607 0 1.000 0.919
glms.lvm.zostera.2 0.154 1.386 1.000 1.662
Parvicardium.exiguum Perinereis.cultrifera
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 1.000 0.824 1.000 0.973 1.000
glms.lvm.zostera.2 0.995 0.063 1.000 2.035 0.992
Perioculodes.longimanus Phoronida Phyllodoce.sp.
Dev Pr(>Dev) Dev Pr(>Dev) Dev
glms.lvm.zostera.1
glms.lvm.zostera.3 2.111 0.995 1.284 1.000 0
glms.lvm.zostera.2 0 1.000 0.11 1.000 0
Platyhelminthes Platynereis.dumerilii
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 1.000 1.372 1.000 1.257 1.000
glms.lvm.zostera.2 1.000 2.602 0.984 0.374 1.000
Polititapes.aureus Polychaeta.larvae
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 0.751 1.000 0 1.000
glms.lvm.zostera.2 6.933 0.404 0 1.000
Polydora.ciliata Polygordius.neapolitanus
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 13.8 0.012 0 1.000
glms.lvm.zostera.2 12.416 0.045 0 1.000
Prionospio.cirrifera Protodorvillea.kefersteini
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 28.158 0.001 7.725 0.294
glms.lvm.zostera.2 4.922 0.749 4.879 0.749
Pseudocuma.longicorne Rissoa.membranacea
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 0.236 1.000 1.543 0.998
glms.lvm.zostera.2 1.386 1.000 1.019 1.000
Rissoa.splendida Salvatoria.clavata
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 5.569 0.607 15.955 0.008
glms.lvm.zostera.2 0 1.000 0 1.000
Schistomeringos.rudolphi Sphaerosyllis.hystrix
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 2.389 0.982 0.023 1.000
glms.lvm.zostera.2 8.799 0.193 4.024 0.863
Spio.filicornis Spisula.subtruncata
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 6.523 0.435 0 1.000
glms.lvm.zostera.2 0.246 1.000 0 1.000
Stenosoma.capito Syllis.gracilis Syllis.hyalina
Dev Pr(>Dev) Dev Pr(>Dev) Dev
glms.lvm.zostera.1
glms.lvm.zostera.3 2.178 0.994 4.755 0.707 2.507
glms.lvm.zostera.2 6.661 0.472 17.831 0.008 0
Tellina.tenuis Thracia.phaseolina
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 0.979 1.622 0.997 0.811 1.000
glms.lvm.zostera.2 1.000 1.386 1.000 0 1.000
Tricolia.pullus Tritia.neritea Tritia.reticulata
Dev Pr(>Dev) Dev Pr(>Dev) Dev
glms.lvm.zostera.1
glms.lvm.zostera.3 5.332 0.629 0.632 1.000 0
glms.lvm.zostera.2 1.38 1.000 3.479 0.917 8.301
Turbellaria Upogebia.pusilla
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 1.000 1.392 0.999 5.769 0.566
glms.lvm.zostera.2 0.239 1.125 1.000 0 1.000
Arguments:
Test statistics calculated assuming uncorrelated response (for faster computation)
P-value calculated using 999 resampling iterations via PIT-trap resampling (to account for correlation in testing.
Well this is tough to interpret.. Multivariate test table’s Dev is decrement from upper model, so each p-value indicates the difference between the model and upper one is statistically significant… But no info on which model represents the species matrix best.
I’ll go with the model with the lowest AIC, since there is no other objective criterion to go with.. This happens to be model 2 (groups = stations). Exactly the same result as from the classical methods.
glms.lvm.zostera.1$AICsum
[1] 6804.496
glms.lvm.zostera.2$AICsum
[1] 6661.507
glms.lvm.zostera.3$AICsum
[1] 6730.346
Now, let’s try to see a different thing - which environmental parameters best describe the species response.
I’m going to use the PCA-filtered environmental data - it’s still going to be a slog, with 7 potential predictors..
First, construct the formula for the model - will do it separately in case I need to update it later, etc.
NB there is year here - I want to try with it first!! And I don’t want the Secchi depth - it has NAs for Vromos.
(formula.env.glms.zostera <- formula(paste("zoo.mvabnd.zostera ~",
paste(env.zostera %>% select(-c(station, secchi)) %>% names(), collapse = "+")))
)
zoo.mvabnd.zostera ~ year + Ntotal + chl_a + LUSI + TOM + moisture_content +
mean_grain_size + sand + silt_clay + shoot_count + ag_biomass_wet +
bg_biomass_wet
Fit the GLMs - including all environmental parameters - to the zostera abundance data.
env.glms.zostera <- manyglm(formula.env.glms.zostera,
data = env.zostera,
family = "negative.binomial")
Explore the fit (residuals, diagnostic plots, etc.).
## residuals vs fitted values
plot(env.glms.zostera)
## all traditional (g)lm diagnostic plots
plot.manyglm(env.glms.zostera, which = 1:3)
# ### source mvabund GLM plotting functions modified to use a grey palette - I just can't redo these plots on my own, the function is doing too complicated things internally to scale the x and y axes
# source(here(functions.dir, "default.plot.manyglm_grey.R"))
# source(here(functions.dir, "plot.manyglm_grey.R"))
#
# par(mfrow = c(2,2))
# lapply(1:3, function(i) plot.manyglm.grey(glms.lvm.zostera, which = i, sub.caption = ""))
# par(mfrow = c(1, 1))
Well, it’s good enough if you ask me (still the kinda strange “line” at lin.pred = -6; otherwise residuals are random enough).
Save the model!
write_rds(env.glms.zostera,
here(save.dir, "glms_env_zostera.RDS"))
Before going off and running an ANOVA to check which predictors best explain the species abundance patterns, I’ll try to reduce the model a little - it might even improve the fit, not to mention the run time.
top.env.glm.red.zostera <- evaluate_glms_env(env.glms.zostera)
the condition has length > 1 and only the first element will be used
[1] "Best model: zoo.mvabnd.zostera ~ Ntotal + sand + shoot_count + ag_biomass_wet + bg_biomass_wet"
Check its fit.
## residuals vs fitted values
plot(top.env.glm.red.zostera)
## all traditional (g)lm diagnostic plots
plot.manyglm(top.env.glm.red.zostera, which = 1:3)
I think it’s fine; might even be better than the full model.. Save it, too.
write_rds(top.env.glm.red.zostera,
here(save.dir, "glms_top_env_red_zostera.RDS"))
Run ANOVA on this model.
(top.env.glm.red.zostera.aov <- anova.manyglm(top.env.glm.red.zostera,
test = "LR", p.uni = "adjusted",
nBoot = 999, ## limit the number of permutations for a shorter run time
show.time = "all")
)
Resampling begins for test 1.
Resampling run 0 finished. Time elapsed: 0.00 minutes...
Resampling run 100 finished. Time elapsed: 0.13 minutes...
Resampling run 200 finished. Time elapsed: 0.26 minutes...
Resampling run 300 finished. Time elapsed: 0.39 minutes...
Resampling run 400 finished. Time elapsed: 0.53 minutes...
Resampling run 500 finished. Time elapsed: 0.66 minutes...
Resampling run 600 finished. Time elapsed: 0.79 minutes...
Resampling run 700 finished. Time elapsed: 0.92 minutes...
Resampling run 800 finished. Time elapsed: 1.05 minutes...
Resampling run 900 finished. Time elapsed: 1.18 minutes...
Resampling begins for test 2.
Resampling run 0 finished. Time elapsed: 0.00 minutes...
Resampling run 100 finished. Time elapsed: 0.14 minutes...
Resampling run 200 finished. Time elapsed: 0.28 minutes...
Resampling run 300 finished. Time elapsed: 0.43 minutes...
Resampling run 400 finished. Time elapsed: 0.57 minutes...
Resampling run 500 finished. Time elapsed: 0.71 minutes...
Resampling run 600 finished. Time elapsed: 0.87 minutes...
Resampling run 700 finished. Time elapsed: 1.02 minutes...
Resampling run 800 finished. Time elapsed: 1.16 minutes...
Resampling run 900 finished. Time elapsed: 1.31 minutes...
Resampling begins for test 3.
Resampling run 0 finished. Time elapsed: 0.00 minutes...
Resampling run 100 finished. Time elapsed: 0.15 minutes...
Resampling run 200 finished. Time elapsed: 0.29 minutes...
Resampling run 300 finished. Time elapsed: 0.43 minutes...
Resampling run 400 finished. Time elapsed: 0.56 minutes...
Resampling run 500 finished. Time elapsed: 0.70 minutes...
Resampling run 600 finished. Time elapsed: 0.85 minutes...
Resampling run 700 finished. Time elapsed: 0.99 minutes...
Resampling run 800 finished. Time elapsed: 1.13 minutes...
Resampling run 900 finished. Time elapsed: 1.27 minutes...
Resampling begins for test 4.
Resampling run 0 finished. Time elapsed: 0.00 minutes...
Resampling run 100 finished. Time elapsed: 0.14 minutes...
Resampling run 200 finished. Time elapsed: 0.27 minutes...
Resampling run 300 finished. Time elapsed: 0.42 minutes...
Resampling run 400 finished. Time elapsed: 0.55 minutes...
Resampling run 500 finished. Time elapsed: 0.69 minutes...
Resampling run 600 finished. Time elapsed: 0.83 minutes...
Resampling run 700 finished. Time elapsed: 0.96 minutes...
Resampling run 800 finished. Time elapsed: 1.09 minutes...
Resampling run 900 finished. Time elapsed: 1.22 minutes...
Resampling begins for test 5.
Resampling run 0 finished. Time elapsed: 0.00 minutes...
Resampling run 100 finished. Time elapsed: 0.12 minutes...
Resampling run 200 finished. Time elapsed: 0.25 minutes...
Resampling run 300 finished. Time elapsed: 0.37 minutes...
Resampling run 400 finished. Time elapsed: 0.50 minutes...
Resampling run 500 finished. Time elapsed: 0.62 minutes...
Resampling run 600 finished. Time elapsed: 0.74 minutes...
Resampling run 700 finished. Time elapsed: 0.86 minutes...
Resampling run 800 finished. Time elapsed: 0.98 minutes...
Resampling run 900 finished. Time elapsed: 1.11 minutes...
Time elapsed: 0 hr 6 min 45 sec
Analysis of Deviance Table
Model: manyglm(formula = zoo.mvabnd.zostera ~ Ntotal + sand + shoot_count +
Model: ag_biomass_wet + bg_biomass_wet, family = "negative.binomial",
Model: data = env.zostera)
Multivariate test:
Res.Df Df.diff Dev Pr(>Dev)
(Intercept) 31
Ntotal 30 1 299.7 0.001 ***
sand 29 1 272.1 0.001 ***
shoot_count 28 1 355.3 0.001 ***
ag_biomass_wet 27 1 235.2 0.001 ***
bg_biomass_wet 26 1 424.9 0.001 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Univariate Tests:
Abra.alba Abra.sp. Actiniaria Alitta.succinea
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
Ntotal 6.177 0.545 4.157 0.945 1.004 1.000 19.751
sand 19.42 0.004 0 1.000 0.412 1.000 0.204
shoot_count 0.506 1.000 0.002 1.000 2.581 0.999 1.336
Ampelisca.diadema Amphibalanus.improvisus
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
Ntotal 0.003 5.846 0.578 10.874 0.057
sand 1.000 25.503 0.002 1.185 1.000
shoot_count 1.000 27.113 0.001 0.224 1.000
Ampithoe.sp. Anadara.kagoshimensis Apherusa.bispinosa
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
Ntotal 10.408 0.080 1.004 0.999 0.83
sand 1.162 1.000 0.381 1.000 3.329
shoot_count 0.23 1.000 2.765 0.995 0
Apseudopsis.ostroumovi Bittium.reticulatum
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
Ntotal 1.000 5.466 0.634 1.145 0.999
sand 0.981 5.211 0.802 0.291 1.000
shoot_count 1.000 10.604 0.169 18.053 0.006
Brachynotus.sexdentatus Capitella.capitata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
Ntotal 0.004 1.000 1.871 0.998
sand 0.202 1.000 1.595 0.999
shoot_count 0.445 1.000 15.66 0.019
Capitella.minima Chamelea.gallina Chironomidae.larvae
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
Ntotal 1.126 0.999 0.789 1.000 9.648
sand 0.223 1.000 3.252 0.986 0
shoot_count 20.087 0.003 12.142 0.098 0.001
Cumella.limicola Cumella.pygmaea
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
Ntotal 0.162 0.501 1.000 2.1 0.996
sand 1.000 2.908 0.992 7.78 0.433
shoot_count 1.000 3.986 0.934 0.001 1.000
Cytharella.costulata Diogenes.pugilator Eteone.flava
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
Ntotal 0.849 1.000 0.003 1.000 4.157
sand 0.163 1.000 0.317 1.000 0
shoot_count 0.48 1.000 6.108 0.672 0.002
Eunice.vittata Eurydice.dollfusi
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
Ntotal 0.940 3.083 0.990 2.1 0.996
sand 1.000 3.191 0.988 7.638 0.449
shoot_count 1.000 0.665 1.000 0.142 1.000
Exogone.naidina Gastrosaccus.sanctus
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
Ntotal 7.771 0.381 0.83 1.000
sand 3.266 0.986 3.329 0.976
shoot_count 0.698 1.000 0 1.000
Genetyllis.tuberculata Glycera.sp. Glycera.tridactyla
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
Ntotal 1.232 0.999 1.097 0.999 0.066
sand 6.164 0.679 4.521 0.887 2.045
shoot_count 3.149 0.976 0.675 1.000 8.534
Glycera.unicornis Harmothoe.imbricata
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
Ntotal 1.000 1.004 1.000 0.915 1.000
sand 0.999 0.381 1.000 6.699 0.595
shoot_count 0.342 2.765 0.996 10.277 0.192
Harmothoe.reticulata Heteromastus.filiformis Hirudinea
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
Ntotal 2.036 0.997 2.133 0.996 0
sand 8.582 0.320 4.357 0.912 0.091
shoot_count 0.077 1.000 1.732 1.000 0.014
Hydrobia.acuta Hydrobia.sp. Iphinoe.tenella
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
Ntotal 1.000 1.425 0.999 1.301 0.999 0.144
sand 1.000 0.408 1.000 0.215 1.000 3.241
shoot_count 1.000 2.507 0.999 0.013 1.000 0.005
Kellia.suborbicularis Lagis.koreni
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
Ntotal 1.000 2.726 0.990 0 1.000
sand 0.986 1.702 0.999 2.559 0.996
shoot_count 1.000 8.058 0.409 11.219 0.140
Leiochone.leiopygos Lentidium.mediterraneum
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
Ntotal 2.493 0.992 1.659 0.999
sand 1.08 1.000 6.526 0.645
shoot_count 7.538 0.473 0.131 1.000
Lepidochitona.cinerea Loripes.orbiculatus
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
Ntotal 0.83 1.000 1.634 0.999
sand 3.215 0.988 1.236 1.000
shoot_count 0.113 1.000 22.136 0.001
Lucinella.divaricata Magelona.papillicornis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
Ntotal 2.519 0.992 0.83 1.000
sand 0.005 1.000 3.215 0.988
shoot_count 3.884 0.942 0.113 1.000
Maldane.glebifex Melinna.palmata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
Ntotal 1.004 0.999 3.974 0.966
sand 0.412 1.000 1.826 0.999
shoot_count 2.581 0.999 5.31 0.810
Microdeutopus.gryllotalpa Micromaldane.ornithochaeta
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
Ntotal 5.334 0.664 4.195 0.900
sand 1.148 1.000 6.449 0.655
shoot_count 3.03 0.981 0 1.000
Micronephthys.stammeri Microphthalmus.fragilis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
Ntotal 3.828 0.970 0.03 1.000
sand 0.614 1.000 1.905 0.999
shoot_count 1.079 1.000 2.619 0.998
Microphthalmus.sp. Monocorophium.acherusicum
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
Ntotal 1.497 0.999 14.561 0.019
sand 0.028 1.000 0.234 1.000
shoot_count 4.309 0.932 18.824 0.005
Mytilaster.lineatus Mytilus.galloprovincialis Nemertea
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
Ntotal 15.75 0.012 1.004 0.999 0.617
sand 1.073 1.000 0.381 1.000 0.254
shoot_count 5.489 0.792 2.765 0.995 9.1
Nephtys.cirrosa Nephtys.kersivalensis
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
Ntotal 1.000 0.458 1.000 1.004 1.000
sand 1.000 0.059 1.000 0.381 1.000
shoot_count 0.288 1.444 1.000 2.765 0.996
Nereis.perivisceralis Nereis.pulsatoria
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
Ntotal 0.83 1.000 2.909 0.990
sand 3.329 0.976 0.077 1.000
shoot_count 0 1.000 0.002 1.000
Nototropis.guttatus Oligochaeta Paradoneis.harpagonea
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
Ntotal 4.526 0.882 5.007 0.722 1.004
sand 1.555 0.999 1.836 0.999 0.412
shoot_count 4.117 0.932 12.359 0.089 2.581
Parthenina.interstincta Parvicardium.exiguum
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
Ntotal 1.000 1.009 0.999 2.023 0.997
sand 1.000 0.001 1.000 0.537 1.000
shoot_count 0.999 10.718 0.169 2.433 0.999
Perinereis.cultrifera Perioculodes.longimanus Phoronida
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept) <NA> <NA> <NA> <NA> <NA>
Ntotal 2.475 0.992 0.055 1.000 0.004
sand 1.989 0.999 0.003 1.000 2.218
shoot_count 0.04 1.000 5.602 0.769 4.34
Phyllodoce.sp. Platyhelminthes
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA> <NA>
Ntotal 1.000 4.157 0.940 1.998 0.997
sand 0.998 0 1.000 0 1.000
shoot_count 0.928 0.002 1.000 0.742 1.000
Platynereis.dumerilii Polititapes.aureus
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA>
Ntotal 1.236 0.999 0.241 1.000
sand 0.216 1.000 6.606 0.624
shoot_count 6.46 0.620 1.073 1.000
Polychaeta.larvae Polydora.ciliata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA>
Ntotal 4.157 0.937 25.833 0.001
sand 0 1.000 10.787 0.131
shoot_count 0.002 1.000 4.281 0.932
Polygordius.neapolitanus Prionospio.cirrifera
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA>
Ntotal 4.157 0.937 3.887 0.967
sand 0 1.000 48.848 0.001
shoot_count 0.002 1.000 0.005 1.000
Protodorvillea.kefersteini Pseudocuma.longicorne
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA>
Ntotal 5.159 0.708 0.884 1.000
sand 0.994 1.000 2.017 0.999
shoot_count 0.73 1.000 0.192 1.000
Rissoa.membranacea Rissoa.splendida Salvatoria.clavata
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept) <NA> <NA> <NA> <NA> <NA>
Ntotal 2.146 0.996 0.416 1.000 1.798
sand 3.924 0.953 3.515 0.964 2.994
shoot_count 0.204 1.000 1.17 1.000 0.118
Schistomeringos.rudolphi Sphaerosyllis.hystrix
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA> <NA>
Ntotal 0.999 1.251 0.999 0.109 1.000
sand 0.991 1.657 0.999 0.151 1.000
shoot_count 1.000 9.19 0.284 10.018 0.223
Spio.filicornis Spisula.subtruncata Stenosoma.capito
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept) <NA> <NA> <NA> <NA> <NA>
Ntotal 17.283 0.007 0.014 1.000 1.53
sand 0.075 1.000 0.081 1.000 4.409
shoot_count 5.161 0.847 3.58 0.961 0.02
Syllis.gracilis Syllis.hyalina Tellina.tenuis
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept) <NA> <NA> <NA> <NA> <NA> <NA>
Ntotal 0.999 3.343 0.986 4.551 0.882 1.144
sand 0.901 0.432 1.000 0 1.000 1.096
shoot_count 1.000 1.132 1.000 0.001 1.000 2.189
Thracia.phaseolina Tricolia.pullus
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA> <NA>
Ntotal 0.999 0.83 1.000 3.903 0.967
sand 1.000 3.329 0.981 0.445 1.000
shoot_count 0.999 0 1.000 0.012 1.000
Tritia.neritea Tritia.reticulata Turbellaria
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA> <NA> <NA>
Ntotal 0.972 1.000 0.049 1.000 0.01 1.000
sand 0.698 1.000 5.174 0.808 0.763 1.000
shoot_count 0.444 1.000 0.182 1.000 0.158 1.000
Upogebia.pusilla
Dev Pr(>Dev)
(Intercept) <NA> <NA>
Ntotal 10.023 0.128
sand 0 1.000
shoot_count 0.001 1.000
[ reached getOption("max.print") -- omitted 2 rows ]
Arguments:
Test statistics calculated assuming uncorrelated response (for faster computation)
P-value calculated using 999 resampling iterations via PIT-trap resampling (to account for correlation in testing.
Nice, all terms are significant now! Again, as with the sand samples, there is one eutrophication term + one sediment composition + seagrass parameters. Seems reasonable enough.
Save this ANOVA (takes too long to run anew).
write_rds(top.env.glm.red.zostera.aov,
here(save.dir, "glms_top_env_red_zostera_anova.RDS"))
Get the taxa with the highest contributions to the tested pattern (here - species most affected by changes in water/environmental quality parameters).
## get the top contributing species for the environmental parameter zostera GLMs
(top.sp.glms.env.red.zostera <- top_n_sp_glm(top.env.glm.red.zostera.aov, tot.dev.expl = 0.75)
)
[1] "Total deviance explained: 0.771"
Polydora.ciliata Alitta.succinea Spio.filicornis
25.833258 19.750751 17.282878
Mytilaster.lineatus Monocorophium.acherusicum Amphibalanus.improvisus
15.750298 14.561059 10.874348
Ampithoe.sp. Upogebia.pusilla Chironomidae.larvae
10.407890 10.022548 9.647785
Exogone.naidina Abra.alba Ampelisca.diadema
7.770666 6.177181 5.846268
Apseudopsis.ostroumovi Microdeutopus.gryllotalpa Protodorvillea.kefersteini
5.466183 5.334367 5.159364
Oligochaeta Syllis.hyalina Nototropis.guttatus
5.007334 4.551098 4.525633
Micromaldane.ornithochaeta Polychaeta.larvae Polygordius.neapolitanus
4.195233 4.157056 4.157056
Eteone.flava Phyllodoce.sp. Abra.sp.
4.157056 4.157056 4.157056
Melinna.palmata Tricolia.pullus Prionospio.cirrifera
3.974488 3.902948 3.886516
Micronephthys.stammeri Syllis.gracilis Eunice.vittata
3.827638 3.343254 3.082611
## unfortunately, mvabund likes to rename my species when converting the data to matrix (no spaces in names), and since I'm going to look them up in my initial untransformed count data, I have to change them back.. DON'T BE IN A HURRY TO DO THAT IF YOU WANT TO SUBSET THE ORIGINAL MATRIX BEFORE RUNNING TRAITGLM
names(top.sp.glms.env.red.zostera) <- names(top.sp.glms.env.red.zostera) %>%
str_replace(pattern = "\\.", replacement = " ")
I’m going to plot these top contributing species, but I’m not using the plot. At least this time it’s more manageable, but still not presentable enough..
## get the species and their abundances from the original count data, and transform them to long format
(abnd.top.sp.glms.env.red.zostera <- zoo.abnd.zostera %>%
select(station, names(top.sp.glms.env.red.zostera)) %>%
gather(key = "species", value = "count", -station) %>%
## turn species into a factor, or you'll be very very sorry later, when they're out of order on the plot. NB need to be in REVERSE order, because ggplot plots from bottom to top, and I want the top-contributing species on top.
mutate(species = factor(species, levels = rev(names(top.sp.glms.env.red.zostera)))) %>%
## add clusters from LVM as a column
mutate(group = case_when(station == "Poda" ~ 1,
station == "Otmanli" ~ 2,
station == "Vromos" ~ 3,
station == "Gradina" ~ 4,
station == "Ropotamo" ~ 5))
)
(plot.top.sp.glms.env.red.zostera <- plot_top_n(abnd.top.sp.glms.env.red.zostera,
mapping = aes(x = species, y = log_y_min(count), colour = factor(group)),
labs.legend = unique(abnd.top.sp.glms.env.red.zostera$group),
lab.y = "Abundance (log(y/min + 1))",
palette = "Set2"
) +
theme(legend.position = "top")
)
Extract the taxon information (univariate tests) from the model ANOVA to present as a table (probably better than this plot, although it’s informative).
## extract the univariate test coefficients (LR) from the environmental model ANOVA. NB keep the row names when converting the matrix to tibble!
table.top.sp.glms.env.red.zostera <- as_tibble(top.env.glm.red.zostera.aov$uni.test, rownames = "var")
## fix the species names - remove first dor
names(table.top.sp.glms.env.red.zostera) <- names(table.top.sp.glms.env.red.zostera) %>%
str_replace(pattern = "\\.", replacement = " ")
## subset only the top species (explaining ~75% of the dataset variation)
table.top.sp.glms.env.red.zostera <- table.top.sp.glms.env.red.zostera %>%
select(var, names(top.sp.glms.env.red.zostera))
## transpose, because a table with 50 columns is just unreadable
(table.top.sp.glms.env.red.zostera <- table.top.sp.glms.env.red.zostera %>%
gather(key = species, value = value, -var) %>%
spread(key = var, value = value) %>%
## arrange as before (terms in the order they appear in the model, and by descending value of the LR for the first model term - here, PO4). Also get rid of the intercept (it's all-NA anyway).
select(species, Ntotal, sand, shoot_count, ag_biomass_wet, bg_biomass_wet) %>%
arrange(desc(Ntotal))
)
Save this to a file - will have to format it as a nice table by hand, unfortunately.
write_csv(table.top.sp.glms.env.red.zostera,
here(save.dir, "taxa_contrib_glms_top_env_red_zostera.csv"))
I’ll do the pseudo-traits analysis - fit single predictive model for all species at all sites, but w/o attempting to explain the different responses using traits - the species ID is used in place of a traits matrix), although I don’t think it will amount to anything useful.
NB only use the top species that exhibited a reaction in the environmental model fit (= the ones accounting for ~75% of the total variability), and only the significant predictors - to improve run times.
sp.response.glms.env.red.zostera <- traitglm(L = mvabund(zoo.abnd.flt.zostera[, names(top.sp.glms.env.red.zostera)]),
R = as.matrix(env.zostera %>% select(Ntotal, sand, shoot_count, ag_biomass_wet, bg_biomass_wet)),
method = "manyglm")
No traits matrix entered, so will fit SDMs with different env response for each spp
sp.response.glms.env.red.zostera$fourth.corner
Ntotal sand shoot_count ag_biomass_wet
names.L.Abra.sp. -0.47066762 0.8734886154 0.46922359 -0.19993821
names.L.Alitta.succinea 0.97052617 0.9239972975 0.15761977 0.34540391
names.L.Ampelisca.diadema 0.07926714 0.4982502059 0.08427469 0.04139468
names.L.Amphibalanus.improvisus 0.01579908 0.1367386026 0.02102261 -0.02672181
names.L.Ampithoe.sp. -0.23185028 0.9766780527 -0.86628249 -0.74571356
names.L.Apseudopsis.ostroumovi 1.64548321 -1.0408187107 -1.81840977 0.56715226
names.L.Chironomidae.larvae 1.31227261 0.1002162364 0.12598021 0.05010550
names.L.Eteone.flava -0.47066762 0.8734886147 0.46922359 -0.19993821
names.L.Eunice.vittata -0.80681453 -1.1651097143 0.87652213 0.87742696
names.L.Exogone.naidina -2.71929348 -1.2204528999 0.50458742 0.72444467
names.L.Melinna.palmata -0.60638932 0.1355595122 0.47942801 -0.75972197
names.L.Microdeutopus.gryllotalpa -0.06793176 0.1548066634 0.05529672 0.15492252
names.L.Micromaldane.ornithochaeta -0.03794287 1.1454709967 0.21794177 -0.18913046
names.L.Micronephthys.stammeri -0.27494238 -0.3866591307 0.52525177 1.09138609
names.L.Monocorophium.acherusicum -0.17196470 0.2856739339 0.13559306 -0.02596581
names.L.Mytilaster.lineatus -0.12863711 0.0415408498 -0.21351554 0.13096055
names.L.Nototropis.guttatus -3.39299319 -2.6671597928 1.95831019 2.81907411
names.L.Oligochaeta 0.24580986 0.0008031048 -0.27535910 0.05999869
names.L.Phyllodoce.sp. -0.47066762 0.8734886150 0.46922359 -0.19993821
names.L.Polychaeta.larvae -0.47066762 0.8734886136 0.46922359 -0.19993820
names.L.Polydora.ciliata 0.08423213 0.0097291405 0.15937837 0.02005751
names.L.Polygordius.neapolitanus -0.47066762 0.8734886143 0.46922359 -0.19993821
names.L.Prionospio.cirrifera 0.88939049 -0.8247531556 -0.04010420 -0.33490099
names.L.Protodorvillea.kefersteini -2.59810569 -2.3545060528 0.75998659 2.29012099
names.L.Spio.filicornis -0.38562919 0.2369199100 0.06576463 -0.07758304
names.L.Syllis.gracilis 0.70617192 -1.4267780844 -1.03624608 0.23992203
names.L.Syllis.hyalina 1.41675539 0.0868009271 0.13451699 0.04413407
names.L.Tricolia.pullus -0.66341187 -2.5031380589 -0.53076627 1.42443911
names.L.Upogebia.pusilla 1.37595789 0.0920536483 0.13121625 0.04647257
bg_biomass_wet
names.L.Abra.sp. -0.25357886
names.L.Alitta.succinea -0.76895635
names.L.Ampelisca.diadema 0.18259939
names.L.Amphibalanus.improvisus -0.23201790
names.L.Ampithoe.sp. 0.03055559
names.L.Apseudopsis.ostroumovi 1.74812763
names.L.Chironomidae.larvae 1.03836845
names.L.Eteone.flava -0.25357886
names.L.Eunice.vittata 0.52390238
names.L.Exogone.naidina 0.09265243
names.L.Melinna.palmata -1.53484185
names.L.Microdeutopus.gryllotalpa 0.19540154
names.L.Micromaldane.ornithochaeta 0.50432844
names.L.Micronephthys.stammeri 0.19869790
names.L.Monocorophium.acherusicum 0.03209065
names.L.Mytilaster.lineatus 0.41368623
names.L.Nototropis.guttatus 1.74403783
names.L.Oligochaeta 0.37862412
names.L.Phyllodoce.sp. -0.25357886
names.L.Polychaeta.larvae -0.25357886
names.L.Polydora.ciliata -0.18809997
names.L.Polygordius.neapolitanus -0.25357886
names.L.Prionospio.cirrifera -0.11604329
names.L.Protodorvillea.kefersteini 1.48926441
names.L.Spio.filicornis -0.32857846
names.L.Syllis.gracilis 1.78601308
names.L.Syllis.hyalina 1.10641092
names.L.Tricolia.pullus 2.08730622
names.L.Upogebia.pusilla 1.07984888
# plot this
a.z <- max(abs(sp.response.glms.env.red.zostera$fourth.corner))
colort <- colorRampPalette(c("blue","white","red"))
plot.spp.z <- lattice::levelplot(t(as.matrix(sp.response.glms.env.red.zostera$fourth.corner)), xlab = "Environmental Variables",
ylab = "Species", col.regions = colort(100), at = seq(-a.z, a.z, length = 100),
scales = list(x = list(rot = 45)))
print(plot.spp.z)
Here at least the directions are a little more coherent than the environmental parameters for the sand stations (seagrass biomasses more or less in the same direction, etc.). The below-ground biomass exerts more pronounced influence on the specific abundances - normal, since most of these are infauna.